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EXECUTIVE  SUMMARY 
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external  heat  flux  from  an  accidental  spill  fire.  This  analysis  considers 
the  assumed  condition  in  which  a  tank  is  completely  engulfed  by  the  flame  of 
a  spill  fire,  the  "worst-case"  condition  resulting  in  a  Boiling  Liquid 
Expanding  Vapor  Explosion  (BLEVE).  Common  examples  for  the  "worst-case" 
situation  include  a  derailed  tank  car  lying  in  a  pool  fire  of  flammable 
material  released  from  a  neighboring  punctured  tank,  or  an  armed, 
mission-ready  aircraft  in  a  shelter  or  on  the  deck  of  an  aircraft  carrier 
exposed  to  an  accidental  jet  fuel  spill  fire. 

This  model  includes  some  of  the  inherent  heat  transfer  phenomena  that 
were  ignored  before.  It  considers  the  thermal  effects  caused  by: 

-  Heat  fluxes  from  total  engulfment, 

-  Initial  transient  heating  of  the  liquid  inside  the  tank  and  formation 
of  a  convective  turbulent  boundary  layer, 

-  Heat  and  mass  transfer  across  the  liquid/gas  interface,  and 

-  High  intensity  heat  input. 

Inclusion  of  the  above  parameters  in  the  heat  transfer  model  has 
improved  the  quality  of  prediction  both  in  the  liquid  phase  and  the  gaseous 
phase.  It  is  recommended  that  efforts  be  continued  to  investigate  other 
geometrical  and  structural  aspects  nf  the  thermal  tankage  problem 
i  nr.  1  ud  i  no  : 

-  Tank  initial  filling  density, 

-  Mechanical  and  physical  integrity  of  the  tank,  and 

-  Structural  response  of  the  tank  shell  to  fire  impingement. 

This  analysis  investigates  and  models  the  response  of  the  tank  under 
thermal  loading  for  two  possible  vented  and  unvented  tank  configurations. 

It  considers  the  thermofluid  dynamic  processes  occuring  inside  the  tank, 
including  the  initial  turbulent  boundary  layer  flow  near  the  wetted  wall  of 
the  tank;  stratification  of  the  tank  bulk  at  the  core;  transport  processes 
at  the  liquid/gas  interface;  convective  heat  transport  in  the  ullage  volume, 
and  transfer  of  heat  fluxes  from  the  bottom  of  the  tank.  Each  of  these 
processes  is  defined,  and  the  governing  equations  are  derived  and  combined 
to  provide  an  integrated  physical  response  model.  The  system  of  governing 
i ntegro/di f ferenti al  equations  representing  dynamics  of  the  physical 
processes  are  solved  numerically,  and  a  model  predicting  the  pressure  rise 
inside  the  tank  is  developed.  The  prediction  code  provides  the  -esponse  of 
a  tank  containing  flammables  t.o  a  large  turbulent  spill  fire  and  estimates 
the  elapsed  time  before  reaching  the  critical  pressure  for  a  given  tank 
configuration.  Also  presented  in  this  report  is  an  evaluation  of  the  code 
performance  and  comparison  of  its  prediction  with  experimental  results. 
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SECTION  I 
INTRODUCTION 


A.  OBJECTIVE 

The  objective  of  this  task  encompasses  two  independent  but  interrelated 
parts:  (1)  development  of  a  model  for  a  large  hydrocarbon  pool  fire,*  and 
(2)  development  of  a  response  prediction  model  for  a  fire-exposed  tank  con¬ 
taining  flammables.  Specifically,  we  are  to  determine  the  conditions  under 
which  a  fire-exposed  tank  containing  flammables  fails  and  Boiling  Liquid 
Expanding  Vapor  Explosion  (BLEVE)  occurs,  and  to  determine  the  elapsed  time 
before  failure.  The  response  of  various  thermophysical  processes  that  occur 
inside  a  tank  as  a  result  of  incident  heat  are  combined  and  an  integrated 
comprehensive  computer  prediction  methodology  is  constructed.  The 
prediction  is  to  provide  the  pressure  rise  history  profiles  for  a  large 
range  of  flammables  and  fire/tank  configurations.  The  prediction  is 
essential  to  on-the-scene  assessment  of  the  fire  environment  and  the  safety 
of  fire-fighting  crews. 

B.  SCOPE/APPROACH 

The  scope  of  this  effort  is  to  develop  an  accurate  heat  transfer  model 
for  the  liquid  and  gaseous  phases  of  a  fire-exposed,  partially  filled  con¬ 
tainer.  In  this  investigation,  the  governing  equations  describing  the  ther¬ 
mophysical  processes  that  occur  inside  a  tank  are  developed  and  solved 
numerically.  Considered  in  this  investigation  are  some  of  the  heat  transfer 
parameters  that  were  ignored  before.  These  include  heat  input  from  the  tank 
bottom  wall  and  initial  transient  heating  of  the  liquid.  The  scope  of  this 
effort  also  requires  conducting  specially  designed  fire  tank  lasts  to  vali¬ 
date  and  assess  the  model  performance. 


*  Development  of  a  large  hydrocarbon  pool  fire  model  has  been  completed.  It 
was  reported  in  a  1986  Task  Report  (NMERI  WA3-12)  entitled  One- 
Dimensional  Hydrocarbon  Fuel  Fire  Model  prepared  by  New  Mexico 
Engineering  Research  Institute  for  AFESC  Engineering  and  Services 
Laborato ry. 
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C.  BACKGROUND 


Large  quantities  of  petroleum  and  volatile  hazardous  chemicals  are 
routinely  extracted,  produced  and  transported,  either  by  hiyhway  or  rail. 

The  possiblity  of  a  major  accident  caused  by  human  error  or  otherwise  during 
transportation  or  maintenance  of  these  facilities  is  quite  real.  The 
impetus  behind  this  analysis  is  to  formulate  and  model  the  thermodynamic 
phenomena  that  occur  in  a  fire-engulfed  tank  and  to  predict  the  elapsed  time 
before  failure  or  the  occurrence  of  8LEVE. 

Before  engaging  in  detailed  mathematics  regarding  modeling  response  of 
a  tank  laden  with  flammables  to  a  1 arge  sooty  fire,  it  is  instructive  to  use 
an  illustrative  example  to  define  trie  sequence  of  events  occurring  in  a 
fi^e/tank  scenario  and  to  demonstrate  the  processes  that  ultimately  lead  to 
3LEVE.  Not  all  fire-engulfed  flammable  tanks  conclude  in  BLEVE.  Consider  a 
large  tank,  e.g.,  a  rail  tank-car  (Figure  1),  partially  filled  with  liquid 
propane  and  the  remaining  volume  (ullage)  occupied  by  propane  vapor  and  air. 
The  ullage  volume  in  the  tank  is  necessary  to  allow  for  swelling  and  expan¬ 
sion  of  the  liquid  that  results  from  incident  external  heat  fluxes. 

When  a  tank  partially  filled  with  volatile  liquid  is  exposed  to  a  large 
turbulent  sooty  fire,  the  heat  flux  from  the  fire  is  partially  stored  in  the 
walls  of  the  tank,  and  is  partially  transferred  to  the  contents,  1  iquid  and 
gas  phases,  of  the  tank.  The  total  heat  transferred  to  the  tank  consists  of 
radiation  from  the  plume  and  convection  from  the  combustion  gases.  The 
contents  of  the  tank  are  heated  by  natural  convection.  Ihe  tank  wall  next 
to  the  gas  phase  always  has  a  higher  temperature  than  the  walls  wetted  by 
the  liquid  phase.  Therefore,  the  temperature  gradient  sets  up  a  conduction 
heat  transport  process  along  the  wall  of  the  tank. 

The  prevailing  heat  transfer  from  the  tank  wetted  wall  to  the  liquid  is 
by  convection.  The  heat  transfer  from  the  wall  establishes  a  turbulent 
convective  thermal  and  vel oci ty  boundary  layer  in  the  liquid  next  to  the 
wall.  The  superheated  turbulent  boundary  layer  next  to  the  wall  raises  the 
heated  liquid  along  the  wall  and  discharges  into  the  subcooled  main  core  at 
the  liquid  gas  interface,  creating  a  convective  circulatory  motion.  This 
process  continues  until  a  near-uniform  temperature  field  resulting  from  a 
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high  rate  of  mixing  is  established.  At  this  tine,  the  heat  transfer  rate  is 
quite  high,  and  almost  all  of  the  incident  heat  transported  to  the  wall  is 
transferred  to  the  liquid;  very  little  goes  into  incr  sing  the  temperature 
of  the  wall,  and  nucleate  boiling  occurs.  As  a  result  of  this,  the  inner 
wall  temperature  in  the  liquid  wetted  region  rarely  exceeds  the  temperature 
of  the  liquid  by  more  than  10  °C  (Reference  1>.  Also,  there  is  a  transi¬ 
tion  region  between  the  liquid  wetted  wall  and  the  vapor  wetted  wall  which 
will  create  a  two -phase  flow  at  the  gas/liquid  interface.  The  thickness  of 
this  region  has  not  been  determined  and  will  depend  on  the  strength  of  the 
boiling  action.  In  this  region,  the  magnitude  of  the  heat  transfer  rates 
will  be  somewhere  between  those  for  the  liquid  wetted  and  vapor  wetted  sur¬ 
faces,  and  as  a  result  the  wall  temperature  in  this  region  will  also  be 
between  those  for  the  liquid  and  vapor  wetted  wall  regions  (Reference  1). 

Tne  prevailing  heat  transfers  in  the  gas  phase  in  the  ullage  volume  are 
by  convection  and  thermal  radiation.  The  initial  convective  heat  transfer 
rate  in  the  vapor  phase  will  be  low  because  of  the  combined  effect  of  the 
low  thermal  conductivity  of  the  gas/vapor  and  the  low  flow  velocities. 

Initial  thermo!  radiation  will  also  be  small  because  of  the  low  wall  temper¬ 
atures.  As  the  temperature  of  the  walls  in  the  ullage  volume  increases,  the 
energy  transport  from  thermal  radiation  will  also  increase.  The  increase  in 
the  energy  transport  in  the  ullage  volume  is  predominantly  due  to  the  radi¬ 
ation  which  is  proportional  to  the  fourth  power  of  the  wall  temperature. 

The  tank  internal  pressure  rises  as  a  result  of  an  increase  in  the  gas-phase 
temperature  and  surface  vaporization  of  the  liquid.  When  the  temperature  of 
the  wetted  walls  reaches  the  saturation  temperature,  correspond i ng  to  the 
instantaneous  tank  pressure,  vapor  is  generated  rapidly  and  the  pressure 
builds  up  quickly.  At  this  time,  the  liquid  portion  of  the  tank  expands  as 
a  result  of  thermal  expansion.  If  the  tank  becomes  completely  filled  with 
liquid  as  result  of  overexpansion,  rupture  will  occur  when  the  critical 
pressure  is  reached.  Tne  critical  pressure  is  related  to  the  decreasing 
strength  of  the  tank  shell  as  its  wall  temperature  increases. 

An  operational  tank  is  equipped  with  a  pressure  relief  valve  which 
would  open  to  allow  removal  of  mass  from  the  tank,  thus, keeping  the  internal 
pressure  below  the  tank  rupture  pressure.  However,  the  released  vapor  from 
the  relief  valve  mixes  rapidly  with  the  ambient  air,  is  ignited,  and  becomes 
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an  additional  source  of  heat  flux  incident  on  the  tank.  The  additional  heat 
flux  that  impinges  directly  on  the  unwetted  wall,  which  is  now  increased  as 
a  result  of  the  drop  in  the  liquid  level,  deteriorates  the  mechanical 
strength  of  the  tank.  Under  these  conditions,  depending  upon  the  integrity 
of  the  physical  structure  of  the  tank,  the  tank  will  either  safely  continue 
to  vent  the  lading  or  it  will  rupture.  If  the  tank  ruptures,  the  remaining 
liquid  in  the  tank  experiences  a  sudden  depressurization.  Since  the  remain¬ 
ing  liquid  cannot  vaporize  quickly  enough  to  establish  a  stable  saturated 
condition,  the  liquid  becomes  superheated.  If  the  liquid  is  superheated 
sufficiently,  a  BLEVE  results. 

Some  of  the  processes  leading  to  the  BLEVE  have  been  defined  by  various 
researchers.  Reference  2  describes  this  phenomenon  in  the  following  form. 

If  a  saturated  liquid  at  high  pressure  suddenly  experiences  a  rapid  pressure 
decrease,  as  ir,  the  case  of  a  rupturing  tank  car,  then  the  liquid  should 
boil  until  the  temperature  of  the  remaining  liquid  is  reduced  to  the  satura¬ 
tion  temperature  for  the  next  pressure  conditions.  This  requires  sufficient 

nucleation  sites  for  the  boiling  to  take  place,  and  these  sites  do  not  exist 
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the  liquid  to  become  superheated  for  a  brief  period  after  the  depressuri za- 
tion  occurs.  If  the  degree  of  superheating  is  sufficiently  large,  homogen¬ 
eous  nucleation  takes  place  throughout  the  bulk  of  the  liquid  and  the  entire 
liquid  mass  vaporizes  almost  instantaneously.  This  rapid  vaporization 
causes  a  vapor  explosion  which  may  cause  the  tank  to  be  ripped  open  and 
sometimes  propels  tank  fragments  for  large  distances.  In  addition  to  the 
shock  loads  from  the  vapor  explosion,  the  vapor  rapidly  mixes  with  the  sur¬ 
rounding  air  and  ignites,  causing  a  fire  ball  which  loads  nearby  structures 
with  intense  thermal  radiation,  Reid  (Reference  3)  suggests  that  a  BLEVE  is 
the  result  of  homogeneous  nucleation  of  a  metastable  superheated  liquid. 

Consistent  with  the  definitions,  occurence  of  3LEVE,  although  not 
totally  independent  of  a  pressure  relief  valve,  is  mostly  controlled  by  the 
initial  liquid  filling  density  and  by  the  magnitude  of  the  incident  heat 
fluxes.  Various  modes  of  response  for  tanks  containing  flammables  have  been 
studied  by  various  researchers.  A  brief  literature  review  is  presented  in 
the  following. 


D.  LITERATURE  REVIEW 

The  response  of  a  tank  engulfed  in  a  pool  spill  fire  has  been  studied 
under  special  conditions  by  various  researchers.  Roberts  et  al .  (Refer¬ 
ence  4)  conducted  experimental  studies  of  insulated  and  uninsulated 
500-liter  tanks.  The  tanks  were  partially  filled  with  liquid  and  exposed  to 
a  kerosene  pool  fire.  In  each  test  they  measured  the  heat  transfer  rate  to 
the  total  system,  the  contents  of  the  tank,  the  boiling  regime,  and  the  tank 
wall  temperature.  They  demonstrated  that  a  response  time  of  3  minutes  (the 
time  needed  for  the  pressure  relief  valve  to  activate)  for  an  uninsulated 
propane  tank  could  be  extended  to  ranges  of  12  to  90  minutes  for  insulated 
tanks  of  similar  material. 

Venart  et  al .  (Reference  5)  conducted  a  study  of  the  physical  behavior 
of  tanks  containing  liquified  fuel  under  accident  conditions.  They  reached 
a  number  of  conclusions  that  are  believed  to  be  true  only  for  low-input  heat 
flux,  i.e.,  heat  fluxes  under  7  w/cm2.  Among  the  conclusions  reached  were: 
(1)  prior  to  pressure  relief  valve  action,  the  primary  source  of  vapor  gen¬ 
eration  originates  in  the  liquid  boundary  layer,  and  (2)  simple  thermo¬ 
dynamic  and  heat  transfer  models  may  be  used  to  predict  the  pressure 
response  of  such  vessels. 

The  problem  of  evaporation  rates  of  liquified  natural  gas  in  containers 
was  considered  by  Hashemi  and  Wesson  (Reference  6).  They  treat  the  system 
as  a  liquid  layer  heated  from  below  and  develop  an  empirical  relationship 
between  the  boi 1  of f  rate  and  the  supersaturati on*  pressure.  This  equation, 
which  shows  that  the  boil  off  rate  varies  as  the  4/3  power  of  the  super-satu¬ 
ration  pressure,  is  solved  in  conjunction  with  the  energy  equation  for  the 
liquid  in  the  tank,  resulting  in  an  implicit  algebraic  equation  in  terms  of 
evaporation  rate. 


*  The  supersatur  .tion  pressure  of  the  liquid  is  defined  as  the  product  of 
the  average  pressure  at  the  saturation  temperature  and  the  total  difference 
between  t.he  temperature  of  the  bulk  of  the  liquid  and  the  temperature  at  the 
sur  face. 
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Mote  that  the  natural  convection  patterns  for  a  liquid  layer  will  not 
exist  for  elevated  pressure  and  the  model  will  not  yield  realistic  solutions 
for  pressure  rise  cases. 

The  natural  laminar  convection  of  fluids  in  an  enclosed  cylindrical 
geometry  has  been  considered  by  Barakat  and  Clark  (Reference  7).  The 
governing  conservative  equations  pertinent  to  the  fluid  inside  the  tank  were 
simplified  by  the  use  of  stream  functions,  and  numerical  techniques  are  used 
for  the  solution  of  these  equations.  They  demonstrated  that  steep  velocity 
and  thermal  gradients  exist  in  the  vicinity  of  the  sidewall  and,  hence,  a 
boundary  layer  flow  occurs.  Analytical  and  experimental  study  of  transient 
natural  convection  in  a  vertical  cylinder  has  heen  reported  by  Evans  and 
Reid  (Reference  8).  Their  experimental  effort  consisted  of  monitoring 
transient  temperature  and  flow  patterns  inside  an  exposed  tank  partially 
filled  with  a  water-glycerin  mixture.  The  tests  were  conducted  for  a  range 
of  Prandtl  numbers  from  2  to  8,000,  and  for  I/O  ratio  (L  =  mixture  depth, 

D  =  cylinder  distance)  and  Orashof  numbers  from  103  to  101 1  encompassing 
both  laminar  and  turbulent  regimes .  They  considered  the  problem  analyti¬ 
cally.  The  system  was  divided  into  three  regions:  (1)  a  boundary  layer 
region  near  the  wall,  (2)  the  mixing  region  at  the  top  of  the  gas/liquid 
interface,  and  (3)  a  miin  core  region  where  the  bulk  fluid  slowly  falls  down 
in  plug  flow.  The  resulting  equations  are  solved  numerically  to  obtain  the 
temperature  profile  in  enclosed  fluids  subject  to  wall  heating.  It  was 
observed  that  after  measurement  of  the  initial  transient  temperature,  the 
core  temperature  increased  linearly  with  time--a  characteri stic  of  low  heat 
t 1 ux  input. 

Transient  thermal  stratification  in  heated,  partially  filled  horizontal 
cylindrical  tanks  was  considered  by  Aydemire  et  al .  (Reference  Q) .  They 
extended  the  analytical  technique  for  boundary  layer  analysis  to  account  for 
the  presence  of  a  covering  of  aluminum  foil  mesh  on  the  inside  walls  of  the 
tank.  Results  v;ith  and  without  ExploSafe  were  compared  with  experimental 
data  for  Freon-113".  Although  fairly  good  agreement  between  the  experi¬ 
mental  and  analytical  results  is  obtained,  the  procedure  appears  to  be  »alid 
only  for  low-input  heat  fluxes.  Vi rk  and  Venkataramana  (Reference  10) 
studied  the  transient  hoiloff  rates  of  a  liquified  natural  gas  from  large 
industrial  storage  tanks  subjected  to  the  perturbation  of  ambient 
atmospheric  pressure. 
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Their  theoretical  model  closely  tracks  the  procedure  covered  in  References 
6-9.  The  model  demonstrates  that  constant  heat  input  to  the  tank  maintains 
a  steady  boiloff  rate  of  LNG  from  th?  tank.  It  establishes  that  a  natural 
convection  boundary  layer  exists  in  the  vicinity  of  the  wall,  and  that  hot 
fluid,  upon  reaching  the  liquid  surface,  diffuses  radically  inward  toward 
the  center,  losing  part  of  its  mass  and  energy  as  boiloff.  They  also  demon¬ 
strated  tnat  step  changes  in  ambient  pressure  (i.e.,  rise  or  drop)  could 
affect  the  boiloff  rate  by  as  much  as  1000  percent. 

8i rk  and  Oosthuizen  (Reference  11)  outlined  a  computer  model  of  a  rail 
tank  car  with  its  exterior  surface  exposed  to  a  spill  fire.  Their  model  is 
two-dimensional,  capable  of  predicting  tank  pressure,  wal  i  temperatures, 
wail  stresses,  vapor  and  1  ■>  n  u  i  d  mean  temperatures ,  flow  rates  through  the 
relief  valve,  and  liquid  level  as  a  function  of  time. 

belichatsios  (Reference  12)  conducted  analytical  and  experimental 
studies  of  fire-exposed  containers  containing  different  solvents  prior  to 
rupture.  He  modeled  the  thermohyd raul i c  process  occurring  inside  the 
container  in  terms  of  a  set  of  physical  dimensionless  prameters.  He  also 
measured  the  container  internal  pressure  and  temperature  rise  resulting  from 
incident  heat  for  different  solvents  using  55-g  lion  steel  drums.  In  this 
investigation  he  concludes  that,  for  most  solvents,  the  time  for  tank  rup¬ 
ture  is  equal  to  the  time  for  the  tank's  wetted  walls  to  reach  the  solvent 
boiling  temperature. 
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SECTION  II 
MODELING  PRINCIPLES 


The  description  and  sequence  of  events  that  occur  inside  a  fire- 
engulfed  tank  are  illustrated  pictorial ly  in  Figure  2.  Figure  2a  shows  the 
fire  engulfment  of  a  rail  tank-car.  As  shown  in  the  figure,  the  tank  is 
nearly  filled  with  the  flammable  liquid  and  the  remaining  volume  is  occupied 
by  vapor.  The  radiative  and  convective  heat  fluxes  from  the  external  fire 
ar-e  conducted  through  the  wall  into  the  tank  lading.  The  heat  input  raises 
the  temperature  of  the  wall  next  to  the  liquid,  and  a  convective  turbulent 
thermal  and  velocity  boundary  layer  is  established.  With  increasing  time, 
the  temperature  of  the  wall  next  to  the  gas  increases  sufficiently  to 
radiate.  The  radiative  heat  flux,  combined  with  the  convective  current, 
causes  vaporization  of  the  liquid  at  the  interface.  At  this  time,  the 
liquid  expands  and  swells,  further  increasing  the  internal  pressure  of  the 
tank.  If  the  tank  is  equipped  with  a  pressure  relief  valve  and  if  the 
tank's  internal  pressure  exceeds  the  valve  operational  pressure,  the  valve 
will  open  to  allow  the  venting  of  lading.  This  will  reduce  the  internal 
pressure  of  the  tank,  but  at  the  expense  of  increasing  the  surface  of  the 
unwetted  wall  in  the  ullage  volume  (see  Figure  2b). 

Simultaneous  with  the  expulsion  of  lading,  two  important  phenomena  that 
significantly  contribute  to  occurrence  of  the  BLEVE  will  take  place: 

1.  The  vented  combustible  vapor  mixes  with  the  surrounding  air  and 
ignites;  therefore,  the  total  radiative  heat  incident  on  the  tank 
is  increased. 

2.  The  radiation  heat  that  impinges  directly  on  the  continuously 
increasing  unwetted  wall  will  result  in  further  reduction  of  the 
tank  shell  mechanical  strength. 

The  additional  radiative  heat  from  the  unwetted  wall  causes  extensive 
surface  evaporation,  a  phenomenon  which  further  raises  the  pressure  consid¬ 
erably.  At  this  time,  the  wetted  wall  temperature  has  reached  the  satura¬ 
tion  temperature,  and  copious  vapor  from  nucleate  boiling  is  generated. 

With  increasing  time,  the  tank's  internal  pressure  rises  quickly,  and  rapid 
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(a)  External  fire  begins.  Heat  is 
conducted  into  tank  loading. 


(b)  Relief  valve  opens.  Venting  of 
of  lading  increases  the  vapor 
space. 


(c)  Wall  adjacent  to  vapor  space 
heats  up  because  vapor  is  a 
poor  conductor  of  heat. 


(d)  Wall  temperatures  rise  to  a 
point  where  the  steel  weakens 
and  the  wall  ruptures. 


Figure  2.  Pictorial  Representation  of  Bleve  (from  Ref.  11) 


generation  of  vapor  and  further  expulsion  of  lading  occur  (Figure  2c).  At 
this  time,  the  strength  of  the  unwetted  wall  has  deteriorated  to  its 
critical  value  and  the  tank  will  rupture  (Figure  2d).  At  the  moment  of 
rupture  when  the  tank  is  depressurized  to  ambient  pressure,  if  the  remaining 
liquid  in  the  tank  is  at  saturation  temperature,  a  violent  explosion,  BLEVE, 
occurs .* 

The  thermofluid  mechanics  of  the  physical  processes  before  BLEVE  were 
discussed  earlier.  They  are  determined  by  the  following  interacting 
transport  phenomena: 

1.  Heating  of  the  walls  of  the  tank. 

2.  Heating  of  the  gaseous  phase  in  the  ullage  volume. 

3.  Heating  cf  the  liquid  phase. 

4.  Transport  of  heat  and  mass  across  the  liquid/gas  interface. 

The  analysis  of  the  mathematical  model  describing  each  of  these  pro¬ 
cesses  is  presented  in  the  following  sections. 

A.  GOVERNING  TRANSPORT  EQUATIONS  IN  THE  ULLAGE  VOLUME 

Figure  3  depicts  the  cross  section  of  an  arbitrary  tank.  This  figure 
describes  interactive  processes  that  occur  in  the  ullage  volume.  The  para¬ 
meters  indicated  are: 


nrr.o 

^rr,i 
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c ,  z 


q" 


Hc,g 
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reradiation  outward  from  the  tank 

reradiation  inward  from  the  wall 

total  heat  flux  from  the  fire  incident  on  the  tank 

radiative  heat  flux  incident  on  the  liquid 

conductive  heat  flux  incident  on  the  liquid 

convective  heat  transfer  from  the  wall 

total  unwetted  surface  area  in  the  ullage  volume 

liquid  surface  area 


*Fcr  more  explanation  of  the  phenomena  refer  to  the  paragraph  entitled 
''Background"  in  Section  I. 
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:  gas  temperature 

T-  :  liquid  temperature  at  the  interface. 

T  :  unwetted  tank  wall  temperature, 

wg 

In  the  derivation  of  the  governing  equations  describing  the  transport 
process  in  the  gas  phese,  the  following  assumptions  have  been  made: 

1.  The  thickness  of  the  tank  wall  is  very  small  compared  with  its 
other  dimensions.  Therefore,  the  temperature  variation  across  the 
thickness  is  negligible. 

2.  Vapor  and  air  in  the  ullage  volume  behave  as  perfect  gases. 

3.  Temperature  and  concentration  distributions  in  the  ullage  volume 
are  homogenous,  since  the  mixing  occurs  relatively  quickly. 

4.  In  spite  of  the  large  temperature  gradients  in  the  wall  near  the 
liquid  meniscus,  the  total  heat  flow  rate  by  conduction  along  the 
wall  is  assumed  negligible, 

5.  The  heat  flux  imposed  by  the  fire  is  uniformly  applied  on  the  tank 
wal Is. 


B.  UNWETTED  WALL  HEAT  TRANSFER  PROCESS 


The  total  heat  stored  in  the  walls  of  the  ullage  volume  is  given  by  the 
following  thermal  energy  balance  equation: 


dT 


p  C  6 

WWW 


dt 


qf 


Vr  ,o 


’rr  ,i 


(1) 


where  p  ,  C  and  6  are  density,  specific  heat,  and  thickness  of  the  wall, 
w  w  w  K 

respectively.  Tne  reradiation  heat  loss  to  the  outside  is  given  by 
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In  Equation  (2),  T  is  the  ambient  temperature,  a  =  0.1713  x  10"S 

oo 

8TU/ft2hrR4  is  the  Stefan-Bol tzmann  constant,  and  e  =  0.99  is  the 
coefficient  of  emissivity.  Similarly,  the  radiation  heat  losses  to  the 
interior  of  the  tank  impinging  on  the  liquid  surface  reads 
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The  total  convective  heat  transfer  rate  in  the  ullage  volume  is  repre¬ 
sented  by  q"  .  It  can  be  shown  (Reference  13)  that  the  total  convective 
c ,  g 

heat  transfer  for  a  given  surface  coefficient,  i.e.,  h,  is  given  by: 


(4) 


where  h  is  the  tank  surface  heat  transfer  coefficient,  and,  as  defined 

earlier,  T  and  T  are  the  unwetted  tank  wall  temperature  and  gas  tempera- 
wg  g 

ture,  respect ively .  The  surface  heat  transfer  coefficient,  h,  can  be  calcu¬ 
lated  from  the  data  given  in  Reference  13,  and  the  convective  heat  transfer 
equation  becomes 


q"  «  C  (T  -  T  )^3 
c  ,g  g  wg  g' 


(5) 


In  Equation  (5)  the  coefficient,  C  ,  depends  on  the  property  values  of  the 
gas  inside  the  ullage  volume.  The  value  of  the  C  is  defined  in  Section 
IV. 


C.  ENERGY  CUNSERVA1 I  UN  IN  THE  GASEOUS  PHASE 


The  conservation  of  energy  in  the  gaseous  phase  is  a  statement  of 
energy  balance  between  the  internal  energy  of  the  gas  and  the  energy  trans¬ 
port  across  the  boundary  and  the  gas/liquid  interface.  The  conservation  of 
energy  balance  reads 


dE 

dt 


q"  A 
c,g  w 


+  M.h . 

i  l 


(6) 


In  deriving  this  equation,  it  is  assumed  that  gas  inside  tne  ullage  volume 
does  not  absorb  any  thermal  radiation.  The  parameters  in  Equation  (6)  are: 


E  :  is  the  internal  energy  of  the  gases  in  the  ullage  volume 
fo.:  is  the  mass  addition  rate  of  vapor  owing  to  the  evaporation 
from  the  liquid  surface, 
h  :  is  the  enthalpy  of  the  vapor 
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The  energy  equation  in  the  ullage  volune  can  be  written  in  terms  of  air  and 
vapor  mass  fractions.  It  follows  that 

M  -MY 
a  a 

M  =  M  Y  (7) 

v  v 

where 

M  :  mass  of  air  in  the  ul  lage  volume 

M  :  mass  of  vapor  in  the  ullaqe  volume 
v 

Y^:  mass  fraction  of  air  in  ullage  volume 
Y^ :  mass  fraction  of  vapor  in  ullage  volume 

Note  that 

Y  +  Y  =  1  (  ) 

a  v 

and  the  total  mass  of  gases  in  the  ullage  volume  is  given  by 

M  =  M  +  M  total  mass  (9 ) 

a  v 

It  is  evident  that  for  a  closed  tank  configuration 
dM 

<nr ' 0  (10) 

Therefore , 


...  dM 
dM  _  _ v  _  f, 

dt  dt  -  1  i 


(ID 


In  the  open  tank  configuration,  the  total  mass  balance  equation  for  the 
ullage  volume  will  take  the  following  form: 


dM 

oT 


(11a) 
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where  dM/dt  is  the  rate  of  change  of  total  mass  in  the  ullage  volume  and  M 
is  the  rate  of  mass  addition  at  the  interface  due  to  the  boiling  and  vapor¬ 
ization.  The  quantity  M  represents  the  mass  rate  of  outgoing  air  and 
vapor,  collectively.  The  mass  rate  of  change  of  each  species  in  the  ullage 
volume  is  calculated  based  on  the  mass  fraction  of  the  outgoing  gases.  For 
example,  the  mass  rate  of  change  of  air  and  vapor  in  the  ullage  volume  is 
given  by  the  following  equations,  respectively: 


d 

dt 


(Y  M)  = 

a 


-Y  M 


a  0 


(lib) 


M. 

1 


Y  M 
v  o 


(lie) 


Using  Eauations  (11a),  (11b),  and  (11c),  the  rate  cf  change  of  air  mass 
fraction  can  be  calculated  as 


(  1  1  n  \ 
V  / 


It  can  be  shown  from  Bernoulli's  law  that  the  outflowing  mass  is  related  to 
the  pressure  drop  in  the  tank: 


(Ilf) 


where  A  represents  the  vent  cross-sectional  area  and  is  the  coefficient 
of  discharge. 


Considering  the  closed  tank  configuration,  the  total  internal  energy 

term  in  the  gas  phase  can  be  written  as  the  sun  of  the  specific  energy  of 

air,  e  ,  and  vapor,  e  : 
a  v 


E  =  M  e  +  M  e 
a  a  v  v 


(12) 


After  taking  the  time  derivative  of  Equation  (12)  and  substituting  for  the 
mass  fractions  from  Equations  (10)  and  (11),  one  obtains 
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dE 

dt 


de 

M  + 

a  at 


de 

v 

""eft 


+ 


M  .e 
i  v 


(13) 


Substituting  for  from  Equations  (12)  and  (13)  into  Equation  (6)  results 
i  n 


M 

a 


M 


de 

\ 

v  dt” 


(14) 


Because  of  the  high  rate  of  mixing,  the  air/vapor  mixture  in  the  ullage 
volune  remains  at  a  uiiform  and  homogeneous  temperature.  Based  upon  this 
conclusion  and  the  definitions  of  internal  energy,  e,  and  enthalpy  of  ideal 
gas,  h,  one  can  modify  Equation  (14).  We  have 


de 

dT 

de 

dT 

4  _  p 

cTT^  La,v 

_ 9 

dt 

* 

V 

dt 

c  — a 

v  V 

’  dt 

h.  =  C 

T 

1  - 

h  = 

C  T  = 

i  v,p 

1 

V 

v,P  9 

(15) 

(16) 


The  time  rate  of  change  of  gas  temperature,  i.e.,  dWdt,  can  be  obtained 
from  the  energy  equation  (Equation  (14))  after  substituting  for  dea/dt, 
dee/dt  and  h ^  from  Equations  (15)  and  (16). 


It  fol 1 ows  that 


dT  q"  A  -  q"  A  +  M.  [C  (T.  -  TV  )  +  ? ,/  pj 
g  _  Hcg  w _ _  cj.  i  L  v,p^  i  g  v __  v 

dt  M  C  *■  M  C 

a  a,v  v  v,v 

Similarly,  for  the  vented  configuration  the  rate  change  of  gas 
temperature  becomes 


(17) 


dT  / 

-n-9-  =  |MC  +  M  C 
dt  \  a  a , v  v  v,v 


-  M  PV/fM  +  M 
o  \  a  v 


-l 


q  A  -  q",A  +  A.  C  t.  -  C  T 
Mcg  w  Mcl  £  i  v,p  i  v,v  g 


(17a) 
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In  arriving  at  the  above  equation  we  have  made  use  of  the  following 
thermodynamic  relationship 

R  P 

r  r  -  _  =  V 

v,p~  Lv,v  ~  TmoT}^  T 

°v  g 

D.  PRESSURE  HISTORY  IN  THE  ULLAGE  VOLUME 

The  total  pressure  in  the  ullage  volume  at  any  time  is  the  sum  of  the 
partial  pressure  of  air  and  the  partial  pressure  of  the  vapor.  Therefore, 
the  total  pressure  becomes 


P  =  P  +  P 
a  v 

and  the  time  derivative  of  the  total  pressure  becomes 
dP  d?a  dPv 

dT  dT~  +  dT"" 


(18) 


(19) 


Using  iaeai  gas  law  relations,  one  can  write 


M  RT 

p  =  _i_  x 

(mol),  V 

a 


p 


V 


M  RT 

v  v  9 
=  — 


(20a) 


(20b) 


For  a  closed  tank  configuration,  the  ulla^w  volume,  V,  remains  nearly 
constant  and  the  time  rate  of  change  of  total  pressure  yields 


dP 


RT 

_ 

V  ( mo  1 )  v 


RM 

v 

V(mol  )v 


RM 

_ a 

V(mol)  . 

a 


■<f r 


(21) 
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In  the  derivation  of  Equation  (21),  the  identities  describing  a  closed  tank 
conf igurat ion,  Equations  (10)  and  (11)  have  been  used.  For  the  vented 
configuration,  the  internal  pressure  of  the  tank  at  any  time  is  determined 
by 


P(t) 


M  (t)  M  (t) 

-2? - +  - 

(mol)  (Mol) 

a  v 


V 


Here 


(mol)v:  molecular  weight  of  vapor 
(mol),:  molecular  weight  of  air 

a 

R:  universal  gas  constant 


E.  MASS  TRANSFER  ACROSS  THE  INTERFACE 


The  gas  temperature  in  the  ullage  volume  is  much  higher  than  the  liquid 
temperature.  Evaporation  at  the  interface  between  the  liquid  and  the  gas 
occurs  because  of  the  convective  heat  transfer  in  the  gas  and  radiative  heat 
flux  from  the  unwetted  wall  impinging  on  the  liquid  at  the  interface  with 
the  vapor  volume.  By  using  stagnant  film  theory  (Reference  14),  one  can 
find  that  tne  evaporation  rate  per  unit  surface  area  is 


—  Ln 


p.g 


C  (T  -  T.) 

i  +  jmLjl _ l 


-  q"  ./m" 
e  f  f  ^  r , }. 


(22) 


Derivation  of  Equation  (22)  is  shown  in  Appendix  A.  In  Equation  (22),  the 
parameters  are 


hf  :  convective  heat  transfer  coefficient  from  a  horizontal 
surface 

C  :  specific  heat  transfer  of  the  gas  mixture 
p.g 

L  effective  heat  of  gasification 
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p,* 


T. 

1 


eff 


specific  heat  of  liquid 
surface  temperature  of  liquid 
substrate  liquid  temperature 

latent  heat  of  vaporization 

£  +  C  (T.  -  T  ) 

P,r  i  l 


F,  ENERGY  TRANSFER  ACROSS  THE  INTERFACE 


(23) 


The  heat  transfer  across  the  liquid  vapor  interface  is  a  statement  of 
energy  balance  between  the  conduction  and  radiative  heat  transfers  incident 
on  the  surface  and  the  energy  associated  with  the  vaporized  mass.  Assuming 
no  radiation  energy  is  absorbed  by  the  vapor  in  the  ullage  volume,  the 
interface  energy  equation  reads 


'eff 


(24) 


The  radiative  heat  transfer  incident  on  the  liquid  surface,  q"  ,  is 

r ,  £ 

directly  related  to  the  total  inward  radiation  energy  emanating  from  the 
side  walls.  It  follows  that 


q"  .A 
rr,i  w 


where,  as  defined  previously,  q"r,r  .  =  ae(T4^  -  TY) 
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SECTION  III 

ANALYSIS  OF  THERMOHYDRAUI.IC  PROCESSES 
OCCURRING  INSIDE  THE  LIQUID  PHASE 


Heat  from  the  external  fire  is  partially  stored  in  the  wall  next  to  the 
liquid  and  partially  transmitted  to  the  liquid  through  convective  currents. 
The  heat  through  the  vertical  walls  establishes  a  turbulent  natural  convec¬ 
tive  velocity  and  temperature  boundary  layer.  The  analyses  of  physical 
phenomena  considered  in  this  investigation  are: 

1.  Determination  of  boundary  layer  flow  produced  by  natural 
convection. 

2.  Evaluation  of  the  core  temperature  variation. 

3.  Analysis  of  liquid  top  layer  where  the  boundary  layer  turns  hori¬ 
zontally  and  the  flow  descends  into  the  liquid  mass. 


The  kinetic  energy  and  the  momentum  of  the  boundary  1 ayer  near  the 
surface  characterize  the  flow  in  the  core.  They  determine  whether  the  flow 
at  the  tank  core  consists  of  large  mixing  eddies  or  a  slow  stratified  motion 
(Reference  15).  The  thermohydraul ic  processes  that  occur  in  the  liquid 
phase  are  shown  in  Figure  4.  The  parameters  indicated  in  the  figure  are 


6  :  boundary  layer  thickness 

Tc(z)  :  core  temperature  (stratification)  distribution 
T,  u  :  temperature  and  velocity  distribution  in  the  boundary  layer 
(defined  later) 

uniform  downward  counter  flow  velocity 
liquid  height  in  the  tank 
area  and  perimeter  of  the  liquid  surface 


H 

l 

A  .  L 


The  mathematical  models  descriptive  of  the  above  physical  processes  are 
developed  in  the  following  sections. 
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A.  HEATING  OF  THE  WETTED  WALLS 

Analogous  to  the  heating  of  the  wall  in  the  gaseous  volume,  it  is 
assuned  that  the  temperature  across  the  thickness  of  the  wall  is  uniform. 
Then,  the  thermal  balance*  equation  for  the  walls  can  be  written  as 


3T 


pwl'w<Sw 


ws 


3t 


=  q  f  -  q 


rr  os 


cts 


(25) 


In  Equation  (25),  q"C£S  Is  the  convective  heat  loss  to  the  liquid  and  Tws  is 
the  side  wall  temperature. 


Similar  to  the  derivation  of  Equation  (5),  the 

the  re-radiation,  q"  ,  heat  transfer  terms  can  be 
’  rros 


convective,  and 

written  respectively  as 


“  Cs(Ls-  Tc)4/3 
5'Vros  1  (t  -  Tt) 

where  the  coefficient  Cs  depends  on  the  physical 
(Reference  13),  given  as 


(25a) 


properties  of  tne  liquid 


CS  =  pocp  I  vc6  g/^5-3)4  Pr23  1 1/3 

B.  CONVECTIVE  FLOW  FIELD  INSIDE  THE  LIQUID 


(25b) 


A  mathematical  description  of  the  heating  of  the  liquid  by  convective 
currents  is  presented.  This  development  ignores  the  possibility  of  laminar 
boundary  layer  which  may  occur  late  in  time.  It  is  also  assumed  that  the 
core  of  the  vessel  does  not  contain  large  eddies,  and  that  the  boundary 
layer  region  can  be  described  continuously  by  turbulent  boundary  1 ayer 
equations.  The  governing  equations  for  the  boundary  layer  are  based  on  the 
integral  representation  and  are  appropr i ately  coupled  with  the  core  flow. 

In  view  of  initially  weak  buoyant  flow  in  the  boundary  layer,  the  Boussinesq 
approximat i on  has  been  implemented. 

*lt  is  assumed  that  the  flames  surrounding  the  tank  are  radiatively  thin. 
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C.  CONSERVATION  OF  MASS 


The  principle  of  conservation  of  mass  is  applied  to  an  elemental  volume 
of  boundary  layer  (see  Figure  5). 


Z 


Figure  5.  Mass  Balance  For  Boundary  Layer  Elemental  Volume 


The  conservation  of  mass  states  that  the  entrainment  rate  is  equal  to 
the  difference  between  mass  efflux  and  influx  in  the  boundary  layer.  There¬ 
fore  the  mass  rate,  dm,  entrained  in  the  boundary  1 ayer  by  convective  cur¬ 
rents  is  equal  to  the  amount  of  mass  leaving  the  boundary  at  any  cross 
section.  It  follows  that 


6  r  ci  i 

J0  [pudy  +  (pudy)  azJ  -  pudy 


=  dm 


where,  after  simplification,  it  becomes 


(26) 


dm 
d  z 


d  ,6  A 

=  ~  J0  pudy 

dz 


(27 ) 
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D.  BUOYANCY  CONSERVATION 


The  buoyancy  equation  is  a  description  of  the  energy  associated  with 
the  upward  How  in  the  boundary  layer  caused  by  the  temperature  gradient. 

The  integral  form  of  the  buoyancy  equation  is  written  for  an  elemental 
volume  of  boundary  layer  (see  Figure  6).  It  follows  that  tne  rate  of  change 
of  energy  per  unit  length  in  the  boundary  1 ayer  is  equal  to  the  algebraic 
sun  of  the  heat  input  (Reference  16). 


Figure  6.  Energy  Balance  for  Boundary  Layer  Elemental  Volume 


The  governing  buoyancy  equation  for  a  tank  having  an  orthogonal  cross 
section  reads 

q"dz  =  d  (/*  puCpTdyj  -  dm  CpTc  (  !8) 

In  Equation  (28),  dz  is  a  di fferenti al  length  in  the  vertical  direction,  u 
is  tne  velocity  of  the  fluid  in  the  boundary  layer,  and  T  is  the  tempera¬ 
ture  of  the  entrained  core  liquid.  Equation  (28)  can  be  rewritten  as 


Substituting  for  dm/d z  from  conservation  of  mass,  Equation  (27),  one 

obtains 


q"/Cp  "  -gj  (/q  puTdy  j  (  /q  Pudy)-Tc 

Since  the  core  temperature  Tc  varies  with  the  liquid  height,  the  last  term 
in  the  rignt  side  of  Equation  (29)  is  modified  to  reflect  tin's  behavior.  It 
follows  that 

wOo  '  st'otc  ^  -?r /o 

Substitution  of  the  above  relation  into  Equation  (29)  gives  the  buoyancy 
equation 


dT 

Zz  I o  ou(T-Tc)  *  f  -  (Jo  zr 


(30) 


Equation  (30)  can  he  modified  to  a  morp  useful  form  if  it  is  multiplied 
through  by  g/(pTQ).  Here,  g  is  the  gravi tat lonal  acceleration  and  Tq  is 
some  reference  temperature.  This  yields 


Define : 


d  r  6 

"dz  o  ug 


T  -  T 


sdy  ~ 


q"g 


pC  T 
po 


d  /  Tc^  ~To  \  -5  H 

-  9*7  — - Poud^ 


(31) 


T  -Tc(2) 


Ac  = 


T  U)  ~ 

C  0 


a  =  6T0a;  aw  =  6  (Tws  -  Tc)  9 


(31a) 


T  = 

w 


pgq 

pCr 


C£S 


( 5 .3)4Pr2 


1/3 


4/3 


q"  =  q,: 

M  hC£S 
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Also,  after  replacing  the  y  coordinate  with  a  stream  function,  ^ ,  one 
obtai ns , 


»!>  =  p0  Jo  ,jdy  (31b) 

£ 

where  ^  =  p0  jQ  udy  is  the  total  flow  rate. 

Then  the  governing  buoyancy  equation  reads 


3_ 

32 


)  r  -  ij;, 

0  W 


8A,: 

3Z~ 


E.  CONSERVATION  OF  MOMENT  OF  MOMENTUM 


(32) 


Recause  the  wall  friction  parameter  in  the  boundary  layer  region  is  not 
well  known,  especially  for  slow  buoyant  flows,  the  motion  of  the  fluid  in 
the  boundary  layer  is  characterized  by  the  moment  of  momentum  instead  of  the 
conventionally  used  momentum  equation  (Reference  17).  This  method  is  quite 
appropriate  for  flows  attached  to  vertical  surfaces.  Rigorous  derivation 
ana  justification  of  this  method  is  beyond  the  scope  of  this  work.  However, 
observe  that  its  acceptabi 1 ity  may  be  based  on  physical  argument,  intuition, 
and  considerable  recent  research  (References  18,  19,  and  20). 

In  deriving  the  moment  of  momentun  equation,  it  is  assumed  that:(l)  the 
flow  is  a  boundary  layer  type  flow,  (2)  the  average  pressure  in  the  plume  is 
the  same  as  the  ambient  pressure  at  sea  level,  which  is  the  consequence  of 
the  Buossinesq  approximation,  and  (3)  the  flow  is  two-dimensional.  We  start 
with  the  Navier-Stokes  turbulent  boundary  layer  equation  (References  21  and 
22). 


p 


+ 


3_ 

sy 


-p  U  V‘ 


3  U 

u  3  y 


+  ApgeTQ 


(33) 
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In  Equation  (33)  z  is  the  vertical  direction,  y  is  the  distance  normal  to 
the  wall,  u  is  the  viscosity,  and  the  quantity  pu"v"  is  the  Reynolds  Stress. 
The  coefficient  $  is  the  thermal  expansion  coefficient  for  liquid,  and  Tq  is 
a  reference  temperature.  Every  variable  in  Equation  (33)  is  represented  by 
its  Favre  average  (Reference  23): 

^  ,1 
u  =  u  +  u 


U  =  (34) 

P 

where  the  velocity  field  u  is  the  algebraic  sum  of  mean,  u  and  fluctuation 
velocity  u".  The  bar  over  the  symbols  denotes  the  time  or  the  ensemble 
average.  Again  we  make  use  of  stream  function,  i  .e . ,  to  replace  the  y 
coordi nate ; 


P  u  =  p 


o  3  y 


p  v  = 


l 
3  Z 


(35) 


With  this  change  of  variable  [(z,y)  +  (z.'lO]  and  neglecting  gradients  of 
stresses  in  the  direction  of  flow,  the  momentum  equation,  Equation  (33), 
takes  the  form 


3_u 

3z 


i! 

3ip 


+  B  i 


which  is  a  more  useful 


T  =  -P /  PQ  U" 

A  =  Ap  /  p  g 


0  4  1  “ 

form.  In  this  equation 

~'r  +  (w/p)  ^/po)2  “ 


(36) 
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To  obtain  the  moment  of  momentum  equation,  first  integrate  Equation  (36) 
over  41  from  4/  to  4^  (total  flow  rate): 


c'Cv  3U  _  3x 

-ip  3z  'ip  34 


d*  *  5  Vi1 


A  /  Udij, 


(37) 


Since 


u  =  0  @  4)  =  4^ 

T  =  0  @  Ip  ~  Ip 


we  have 


—  r.Z  udip  =  -r  +  sT_  (  z  a  /  ud4 

3  Z  '  \p  v  4.  0  '  41  v 


(38) 


Next  integrate  over  4)  from  0  to  4^  .  One  obtains 


in1  a^T7  /0tud*  "  "  iV  TdlJ/  +  I  o''  A  ! 


t  ,  ~ 


Jo 


(39) 


using  the  Leibnitz  rule  (Reference  24) 


h  C  d*  i'o  ,d*  *  BTo  / 0  d*  lo  4 ' ud» 


which  after  partial  integration  reads 


~  /0t  11^4,  =  ST0  /Qt  a  /  uijjd 4)  -  JQt  Td*  (40) 

Equation  (4C )  is  the  desired  moment  of  momentum  equation.  In  general,  the 
second  term  at  the  right  side  of  Equation  (40)  is  small  relative  to  the 
first  term  at  the  right  side  of  the  same  equation;  its  magnitude  is  slightly 
negative  for  turbulent  flows. 

Equation  (40)  demonstrates  that  the  core  stratification  affects  the 
momentum  of  the  boundary  layer  only  indirectly  through  the  temperature  dif¬ 
ference  A  .  In  Equation  (40),  8  is  the  coefficient  of  thermal  expansion  for 
the  liquid  and  has  the  units  of  inverse  temperature. 
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In  deriving  the  moment  of  momentum  equation,  it  is  assumed  that  the 
core  velocity,  u,  is  significantly  small  indeed,  since  the  width  of  the  tank 
is  much  larger  than  the  boundary  layer  thickness  (see  Figure  4).  Finally, 
the  shear  stress  term  in  Equation  (40),  although  small,  must  be  included  in 
the  present  application  because  the  velocities  and  the  growth  of  the 
boundary  layer  are  slowed  down  by  strati f i cati on .  An  approximate 
development  of  the  shear  stress  is  given  in  Reference  22.  It  reads 


1/2 

;  xdp  =  lOo ( F  v)  0 
0  w  t 


Thus, Equation  (40)  may  he  written  as 


a  -  ~  1/2 

! o  u*d,k  =  /0  4  /  updp  -  10d(Fwv) 


(41) 


where  u  is  the  flow  velocity  in  the  boundary  layer,  n  is  the  kinematic 
viscosity  of  liquid,  and  F^  is  as  defined  in  Equation  (31a). 


F.  WALL  HEAT  FLUX  TEMPERATURE  RELATION 

Analogous  to  the  principles  defined  earlier  for  the  gaseous  phase,  the, 
wetted  wall  heat  transfer  is  based  on  calculation  of  turbulent  free  convec¬ 
tion  over  a  vertical  plate  (Reference  25  and  Reference  13,  pp.  197-207)* 

The  convective  heat  transfer  coefficient,  h  ,  for  a  large  range  of  Prandtl 
numbers  (Pr=  C^u/k)  is  given  by 


Nu  =  C j(Gr.Pr) 1/3 


(42) 


where  the  symbol  Mu  is  the  Nusselt  number,  (Nu  =  h  L/K),  and  Gr  is  tne 
Grashof  number  (Gr  =  g8ATL3p2 /u2 ) .  Using  the  experimental  data  given  in 
Reference  13,  one  obtains  the  following  wall  heat  flux  temperature 
relation: 


A  = 
w 


1/4 


(Pr) 1 /2 


(43) 
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G.  VARIATION  OF  CORE  TEMPERATURE 

The  variation  of  the  core  temperature  has  been  considered  for  heating 
and  nonheating  of  the  tank  bottom  surface.  In  the  absence  of  trie  bottom 
heating,  the  core  temperature  will  increase  because  the  hot  layer  near  the 
interface  moves  aown  the  center  of  the  tank.  In  this  case  the  variation  of 
the  core  temperature  can  be  readily  obtained  from  the  general  diffusion  heat 
transfer  equation.  Neglecting  the  second-order  diffusion  terms  ana  the 
viscous  effects,  the  tank  core  temperature  variation  equation  reads 


3A„ 


3t 


-U(z) 


(44) 


In  Equation  (44)  t  -  pg(Tc~  TQ ) ,  where  3  is  trie  expansion  coefficient  and 
U(z)  is  the  core  main  counter  flow  velocity.  Note  that  U(z)  can  be  obtained 
from  the  following  continuity  argument.  Since  the  boundary  layer  thickness, 
5,  is  much  smaller  than  the  tank  diameter,  the  flow  velocity  in  the  boundary 
layer  must  be  much  larger  than  core  velocity.  Therefore,  based  upon  the 
conservation  of  mass  argument,  the  core  velocity  can  be  calculated  by  set¬ 
ting  the  upward  total  mass  flow  in  the  boundary  eaual  to  the  downward 

'l  £ 

total  mass  flow  -pA^U(z)  in  the  core.  It  follows  that 

U(z)p  Af  -  -  *t(z)  iz  (45) 


If  the  bottom  of  the  tank  is  also  exposed  to  an  external  heat  flux, 
adoitional  heat  is  transferred  to  the  bulk  of  trie  liquid.  In  this  case, 
similar  to  the  previous  analysis,  one  must  consider  the  additional  heat 
transfer  processes: 


1.  Heating  of  the  bottom  "horizontal"  surface  (similar  to  the 
analysis  for  "vertical’1  wall). 


2.  Convective  heat  loss  from  the  bottom  to  the  liquid  (similar  to 
Equat  ion  (5)  .  We  have 


dT 


o  L  6 
w  w  w 


dt 


-•  =  q"f 


q"  .  -  q" 

M  rrob  M  ctb 


(46) 
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where 


"‘cb  =  Cb'rw,b  -  \,^'3 


T  •  wall  temperature  at  the  bottom 

w, b 

T  .  :  core  liquid  temperature  near  the  bottom 

x,  9  D 


The  additional  heat  from  the  bottom  is  distributed  by  the  buoyant 
thermals  near  the  bottom  of  the  tank.  We  assume  that  these  thermals  are 
"buoyantly"  strong  enough  to  transfer  the  heat  uniformly  through  the  bulk  of 
the  fluid  while  the  boundary  layers  remain  unaffected.  Under  this  condition 
the  governing  heat  transfer  equation,  Equation  ('4),  becomes 


-U(z) 


3*c  Ab98 
3z  PCpVT  qc£b 


(47) 


In  Equation  (47),  is  the  bottom  surface  area  and  V-j-  is  the  tank  volume. 


The  following  empirical  expressions  (Reference  25)  describe  the 
profiles  of  velocity  and  temperature  in  the  boundary  layer: 


(48) 


The  governing  equations  that  have  been  developed  for  each  zone  together 
with  the  boundary  layer  velocity  and  temperature  profiles.  Equation  (48), 
constitute  a  set  of  deterministic  i ntegro-di fferenti al  equations  that  are 
solved  ntmerically.  The  solutions  provide  internal  pressure  rise  history 
and  boiling  temperature  which  determines  the  tank  response. 
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SECTION  IV 
NUMERICAL  SOLUTION 


A.  NUMERICAL  SCHEME 

A  special  numerical  scheme  is  developed  to  solve  the  governing  thermo¬ 
hydrodynamic  equations  of  processes  occuring  inside  the  tank.  This  scheme 
utilizes  a  stationary  and  a  moving  grid  that  are  essential  to  temporal  and 
spatial  calculations  of  the  liquid  parameters  under  thermal  load.  The 
numerical  scheme  also  includes  a  rezoning  procedure  necessary  to  model  the 
variation  of  the  temperature  in  the  descending  liquid  inside  the  tank  as  a 
funct  ion  of  height . 

The  governing  partial  differential  Equations  (1),  (10),  (11),  (17), 
(22),  (25),  (32),  (46),  (47)  are  solved,  and  the  height-dependent  variables 
are  approximated  in  either  of  the  two  grid  systems  using  a  time-dependent 
staircase  (i.e.,  piecewise  constant)  function  of  height  z.  The  staircase 
function  is  restricted  to  jump  only  at  the  grid  points.  Figure  7  illus¬ 
trates  an  i i'ib> Idritaueous  plot  of  a  staircase  function  of  tne  core  temperature 
T  (t,z)  in  a  stationary  grid  zq,  z^.-.z^,  as  well  as  a  moving  time-depen¬ 
dent  grid,  z^,  z{(t), . .  .z^(t),  zJ^j.  The  prime  sign  indicates  the  moving 
grid.  The  vertical  axis  measures  the  two  corresponding  staircase  approxima¬ 
tions  for  the  core  temperature  T  (t,z).  The  time  axis,  not  shown  in  this 
figure,  is  perpendicul ar  to  the  plane. 

Equation  (25),  describing  the  temperature  of  the  wetted  wall,  can  be 
easily  solved  by  using  a  stationary  grid.  However,  a  moving  grid  which 
follows  the  motion  of  the  descending  liquid  is  required  to  evaluate  the  core 
temperature  variation.  It  is  important  to  note  that  the  temperature 
T  (t,z),  which  is  a  function  of  time  and  height,  is  the  only  variable  that 
needs  to  be  calculated  in  terms  of  both  grids.  This  is  because  of  the 
convective  heat  transfer  term,  Q"c^s>  in  Equation  (25)  which  depends  on 
T  (t,z).  Therefore,  7'c(t,z)  is  calculated  once  in  terms  of  the  moving  grid 
for  the  core  temperature  equation  and  a  second  time  in  terms  of  the 
stationary  grid  for  the  wetted  wall  equation. 
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A  rezoning  procedure  has  been  devised  to  couple  the  two  grid  systems. 
This  procedure  converts  a  staircase  function  in  a  moving  grid  to  a  staircase 
function  in  a  stationary  grid.  Therefore  the  rezoning  procedure  enables  one 
to  convert  from  one  grid  system  to  the  other.  Details  of  this  procedure  are 
presented  in  the  following  section. 


Tc  (t.z) 

♦ 


z.  zj  (t)  z'z  (t)  zj  (t)  z^  (t)  z'b 

J! _ L_i _ I _ I _ 1 _ |_J _ L 

Z0  Z1  z2  Z3  Z4 


Figure  7.  Staircase  Representation  of  the  Core  Temperature  in 
Terms  of  a  Moving  and  Stationary  Grid,  n  =  4 
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B.  COMPUTATIONAL  GRID 

The  variation  of  the  wetted  wall  temperature  as  a  function  of  height  is 
calculated  numerical ly,  based  on  a  stationary  grid  given  as 


z,  =  kH  /n  for  k  =0, 1, . . .  ,n  . 

k  i 

In  this  equation  the  height  of  liquid,  H  ,  is  divided  into  n  cells  of  equal 

size  in  the  interval  0  <  z  <  H  .  Once  the  stationary  grid  is  set  up,  the 

i  .  .  . 

wetted  wall  temperature  T  can  be  approximated  using  a  staircase  function. 


T 

ws 


( t ,  z ) 


f°r  Vl  <  z  < 
for  k=  1,2,. . .  ,n. 


As  discussed  earlier,  the  staircase  function  is  restricted  to  jumps  only  at 
points  ^  • 

The  variation  of  the  core  temperature,  1  (t,z),  due  to  the  motion  of 
descending  liquid  is  modeled  by  using  the  following  moving  grid, 

0  =  z;  <  z;(t)  <  zj(t)  <  ...<  z;(t)  <  z^+i  =  Ht. 

In  this  expression  the  grid  moves  with  the  liquid,  and  Zg= 0  and  z^+j  =  ar 
the  fixed  end  points  corresponding  to  the  bottom  of  tne  tank  and  top  of  the 
liquid,  respectively.  Then  the  core  temperature  can  be  approximated 
using  a  staircase  function  in  the  moving  grid 


Tjt.z)  =  Tc,k(t)  f0r  Zk-l(t)  <  Z  <  Zk(0 

i or  k=l n+ 1. 

A  similar  approach  is  taken  in  a  related  problem  by  Germeles  (Reference  26). 
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The  wetted  wall  temperature  equations,  Equations  (25)  ana  (25a), 

contain  a  convective  heat  transfer  term  which  depends  on  the  core 

temperature,  T  .  Therefore,  analogous  to  the  discretization  of  the  wetted 

wall  temperature,  T  the  core  is  described  in  the  same  manner  using  a 

ws , 

staircase  function  in  a  stationary  grid.  Hence,  T  in  a  stationary  grid 
becomes 


T 

c 


(t,z) 


”  Tc,k(t> 


f°r  zk_j  <  z  <  zk 
for  k  =  1 ,2, . . . ,n . 


A  rezoning  procedure  is  employed  to  convert  the  staircase  function 
in  the  moving  grid  into  a  staircase  function  approximation  in  a  stationary 
grid.  The  application  of  the  two  grid  systems  for  calculating  t'^-;  liquid 
core  temperature  Tc(z,t)  as  a  function  of  height  is  presented  in  the 
fol lowing. 

Assune  the  following  time  grid 


‘k  =  V  kit 


for  k  =  0,1,2, .. . 


arid  consider  the  time  step  from  t .  to  t 


j+1 


Initially,  at  time  t.,  the  liquid  core  temperature,  T  (z,t;),  is  known  and  is 

J  o 

defined  in  terms  of  a  staircase  function  in  the  stationary  grid.  In  other 

words,  the  values  of  T  ,  (t.)  are  known  for  all  k's. 

c,k  j 

STEP  1.  Define  the  moving  grid  z ^  in  the  same  time  interval,  i.e., 

t  .  <  t  <  t .  , ,  so  that 
J  J+l 

(a)  at  the  start  of  the  time  step,  t  =  t  .,  the  moving  grid  coincides 
with  the  stationary  g- id  for  all  values  of  k,  i.e.. 


(V  =  2k 


z  ..=  z 
n+1  n 


for  k  =  0,1 ,2, ...  ,n 
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(b)  the  moving  grid  points  z[  (t)  ,zi>  (t) , . . . .  zPt)  move  with  the  liquid 
core,  while  the  endpoints  z^  and  z^+1  remain  fixed.  Note  that  the  moving 
grid  is  redefined  at  the  start  of  each  time  step. 

STEP  2.  Next,  express  the  core  temperature  T,(z,t)  at  time  t.  in  terms 

J 

of  a  staircase  function  in  the  moving  grid  .  Since  the  moving  grid  coin¬ 
cides  originally  with  the  stationary  grid,  it  follows  that 

=  Tc,k(tj>  for  k  ■  1'2 . "• 

STEP  3.  Iterate  the  discretized  differential  equations  for  T  (z,t) 
over  time  to  calculate  the  staircase  expression  for  the  core  temperature  in 
the  moving  grid,  i.e.,  T'  .  ( t , )  for  k  -  l,2,...,n+l. 

C ,  K  J  + 1 

STEP  4.  The  discretized  side  wall  temperature,  T^,  Equation  (25),  is 

solved  numerically.  This  iterative  procedure  takes  the  staircase  represen- 

t.atici  of  T  at  time  t.  in  a  stationary  grid  and  calculates  the  expression 
ws  J 

for  T  at  the  next  time.  t..,.  The  values  of  T.  ,  ft.)  are  needed  for  this 
vs  '  ji-j.  C,K  j 

calculation  of  the  side  wall  temperature. 

STEP  5.  Using  the  rezoning  method,  the  staircase  function  for  Tc(z,t) 
at  ti-ie  tj  +  j  in  the  moving  grid  is  converted  into  a  staircase  function  in 
the  stationary  grid.  Therefore,  the  values  of  P  ^(P+j)  are  used  for 
calculating  the  values  of  Tc  ^(t-+^)  for  k  =  l,2,...,n  in  the  stationary 
grid. 


This  procedure  is  repeated  to  calculate  the  variables  for  the  subse¬ 
quent  t ime  steps , 

A  detailed  description  of  the  rezoning  method  is  presented  in  the  fol¬ 
lowing.  Consider  the  kth  stationary  cell  z^_^<  z  <  z^,  and  assinie  the 
simple  case  where  it  overlaps  only  with  kth  moving  cell  z  <z^  and  the 

(k+l)th  moving  cell  z£,  <  z  <  z^.+^.  P  other  words,  one  has  the  situation 
where 


Pi  <  Pi  <  <  2k  <  Pr 
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Then  the  core  temperature,  T  ,  in  tne  kth  stationary  cell  given  by  the 

L  j  K 

rezoning  method  is 

Tc,k  =  <2k  -  Jk-I>'1  [<n  -  Vl>  Ti,k  *  (2k  -  2PTi,k+l] 

In  the  above  relation  T  .is  a  weighted  average  of  the  temperatures  T1 

c  >  K  c  ,  k 

and  T^  k+^  in  the  kth  and  (k+l)th  moving  cells,  respectively.  The  weights 
are  proportional  to  the  lengths  of  overlap  z^  -  zk  1  and  zk  -  zk  of  the 
stationary  cell  with  the  moving  cells.  Figure  8  shows  the  relationship 
between  the  moving  and  stationary  grid  systems. 
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Figure  8:  Moving  and  Stationary  Grid  Systems 
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C.  MODELING  TEMPERATURE  OF  THE  WETTED  WALL 


The  governing  temperature  equation  for  the  side  wall  next  to  the  liquid 

T  (t.  z)  was  derived  earlier  in  Equation  (25).  It  satisfies  the  following 
w  s 

partial  differential  equation 


3t  e 

r  r  ws  _  *1,  *it  *u 

pwSv5w  “at  Q  f  “  Q  rros  -  q  C£S 


(49) 


where 


q"  =  ea  (T  **  -  T  M 
Hrros  '  ws  ®  ’ 


q"  =  C  (T  -  T  )4'3 
Mc  jts  s  v  ws  c 


Cs=  Po  CP 


1/3 


L (5.3)**Pr2  J 

,  -4/3 

Clv=  (5.3)  =  0.11 


=  C 


vcfig  \  1/3 


*v  0  P\  pr 2 


The  aoove  partial  differential  equation  contains  a  derivative  expli¬ 
citly  with  respect  to  time  t  but  not  z.  Therefore,  it  can  be  replaced  by  a 
system  of  ordinary  differential  equations  involving  ordinary  derivatives 
with  respect  to  t.  Furthermore,  the  z  dependent  variables  in  Equation  (49) 
are  all  replaced  by  their  corresponding  staircase  functions  in  a  stationary 
grid. 


The  governing  expressions  for  the  z  dependent  variables  as  staircase 
functions  in  the  stationary  grid  are 


TWS  «-•*>■  TWS,k  l‘> 

for 

zk-l 

Tc  o.  a  ■  Tc>k  <t) 

for 

zk-l 

q"  (t,  z)  =  q"  (t) 

rros  rroSfk 

for 

Zk-1 

q"  (t,  z)  =  q"  .  (t) 

CHS  *  Hczs,k  v  ' 

for 

Zk-1 

39 


and  the  partial  differential  equation  for  T  (t,  z)  in  the  staircase 
function  approximation  becomes  the  following  system  of  ordinary  differential 
equations , 


dT 


p  C  5 

WWW 


dt  q 


rros  ,k 


-  q 


cts  ,k 


(50) 


where 


q"  .  =  eo  (T  4  -  T  4) 
rros,k  v  ws,k  °° 


q"  ,  ,  -  C  (T  ,  -  T  .  ) 

^c£s,k  s  ws,<  c,k 


4/3 


for  k  =  1,2, ...  ,n. 


The  ordinary  differential  equation  for  the  temperature  T  ^(t)  of  the 
bottom  wall,  liquation  (46),  reads 


dT  u 

~  _  wb  •  „  • , 

p  c  5  ,  —j-p--  =  q  *  -  q 

wwwdt  4  f  ^ 


rrob  "  q  ctb 


(51) 


where 


q"  ,  =  eo  (T4.  -  T  4) 

Hrrob  '  wo  *  ' 


Hc£b 


=  Cb  <Twb  -  Tcb^/3 


/V9 

C,  =  Clup  C  -- — 

b  lh  0  p\  Pr2 


1/3 


=  (C ../  C,  )C 
v  lh  lv  s 


C,.  =  0.06 
lh 

The  constants  C,,  and  C,  are  coefficients  within  the  heat  transfer  coeffi- 
lh  lv 

cients  correspondi ng  to  the  horizontal  and  vertical  walls  of  the  tank,  res¬ 
pectively.  The  quantity  T  is  the  liquid  core  temperature  at  the  bottom. 
The  quantity  T  (t)  equals  the  value  of  Tc(t,z)  in  the  first  step  of  its 
staircase  approximation. 


Tcb  <*“>  *  Tc,I  (t> 
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A  fourth-order  Runge  Kutta  nunerical  method  (Reference  27)  is  used  for 
solving  the  set  of  ordinary  differential  equations.  Equations  (50)  and  (51). 
The  quantities  are  kept  const  *nt  during  the  calculation  of  a  single 

time  step.  However,  they  a»'e  up<  ated  to  new  values  at  the  end  of  each 
iteration. 

0.  MODELING  OF  LIQUID  CORE  T  .MPERATURE 


The  partial  differential  equations  for  the  liquid  core  temperature 
T_(t,z)  was  derived  earlier  t  Equation  (47);  it  reads, 


3T  3T  A 

_ S.  +-  u _ M _ q« 

at  u  az  P  cpvT  qc£b 


(52) 


where  V  is  the  tank  total  liquid  volume  and  the  vertical  velocity  U(t,  z)  of 
the  core  is  given  by  Equation  (45)  as, 

U ( t ,  z)  =  -4t(t,  z)  L ^  /  ( pqA  J  (53) 

Here  L  is  the  perimeter  of  the  horizontal  cross  section  and  \K.(t,  z)  is  the 
l  L 

boundary  layer  mass  flux.  Detailed  description  of  \jit(t,z)  is  presented 

in  the  next  section.  The  left  side  of  the  Equation  (52)  is  the  material 

derivative  OT  /Dt  (also  called  total  derivative)  of  T  with  respect  to  time, 
c  c 

Therefore,  for  a  grid  which  moves  with  the  liquid,  the  total  derivative  is 
replaced  in  the  discretized  equations  by  an  ordinary  derivative  dTp  ^/dt. 
Hence,  the  discretized  expression  for  T  (l,  z)  as  a  staircase  function  in 
the  moving  grid  becomes 


Tc(t,  z)  =  T^k  (t)  for  zj_x(t)  <  z  <  zMt) 

for  k  =  1,  2,...,n  +  1. 

and  the  governing  partial  differential  equation,  Equation  (52),  in  the 
staircase  function  approximation  will  become  the  following  system  of 
ordinary  differential  equations. 


dT'  L  Ak 
c,k  -  b  •„ 

~ar~  P  c  \rT  q  cub 


for  k  =  1,2,  . . ,  ,n  (54) 
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Equation  (52)  must  satisfy  the  boundary  condition  at  the  top,  z  = 


T  (t,H  )  =  T 
c'  *  i  c 


top 


( t) 


The  quantity  ,  core  temperature  at  the  top,  is  determined  by  the  rate 

that  the  mass  and  heat  ^"e  being  added  tc  t tie  core  top  from  tne  convective 
boundary  layer  flow.  C  .''responding  to  the  transformation  of  partial  differ 
ential  Equation  (52)  into  tiio  system  of  ordinary  differential  equations 
shown  in  Equation  (54)  by  application  of  the  staircase  approx  imat  ion,  the 
above  boundary  condition  will  take  the  following  form: 


TU i(t)  * 


Tc  too' 0 • 


The  variables  appearing  in  the  liquid  core  velocity  equat ion,  Equation  (53) 
are  defined  in  a  stationary  grid  as: 


U,  ( t )  -II  ( t ,  z.  )  for  k  =  0,1 _ n  . 

ft  K 

Note  that,  unlike  the  core  temper ature ,  $t(t , z)  and  U(t,z)  are  not  being 
represented  as  staircase  functions,  and  are  not  constant  between  the  grid 
points.  Then  Equation  (53)  at  the  grid  points  becomes 

Uk(t)  *  •«'t,k(t,Lt/  <3oV  (55) 

Euler's  method  (Reference  27)  is  used  for  numerically  updating  the 
values  of  T‘  ,  during  each  Lime  increment  At.  Also,  the  values  of  the  mov- 

j  K 

ing  grid  points  at  the  end  of  a  time  increment  for  all  values  of  k  are 
calculated  from: 


z£(tj+  At)  =  Zj^  it.)  ^  At  Uk(tj.)  =  Zk  +  At  lJk(t.) 


and 


n+l  =  V 
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The  rezoning  procedure  is  used  to  convert  T(  ^  in  the  moving  grid  into 
T  i  in  3  stationary  grid  at  the  end  of  each  time  step. 

C  j  K 

E.  MODELING  OE  CONVECTIVE  BOUNDARY  LAYER 

In  the  convective  boundary  layer,  a  continuous  upward  flow  of  mass  and 
heat  takes  place.  The  motion  of  the  mass  in  the  boundary  layer  is  defined 
by  the  buoyancy  equation.  Equation  (29).  The  buoyancy  equation  was  recast 
into  the  more  useful  form  in  terms  of  total  mass  flux  \|i  in  Equation  (32). 

In  a  similar  fashion,  the  upward  heat  flux  is  defined  and  expressed  in  terms 
of  >|)  .  The  values  of  mass  and  heat  flux  at  the  top  of  the  boundary  are  used 
to  estimate  the  temperature  T^  (t)  at  the  top  of  the  core.  The  buoyancy 
equation  together  with  the  temperature  profile  (Reference  28)  and  the  heat 
flux  temperature  relations,  in  the  boundary  layer,  are  given  in  Equations 
(56),  (57),  and  (58),  respectively. 


where 


Z  =  eg  (T  -  T  ) 
w  v  ws  c- 


Integrating  Equation  (57)  from  p  =  0  to  ip  =  ^  gives. 


J0  Ad*  =  0.937  aw  Pr-i'2  vc1/4  <|>t3 


/4 


(59) 


Using  Equations  (58)  and  (59),  the  buoyancy  equation,  Equation  (56),  reduces 
to 

aL 


~  (A  <1,  3/4)  =  C  A  4/3  -  C.  ip.  — 1 
3Z  '  W  yt  '  0  W  1  Yt.  D  Z 


(60) 
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where 


Co  =  P0vc1/12/[(‘937)(5*3)  4/3Prl/6] 

C,  =  Prl/2/^0.937)  vc1/4 

This  is  a  form  of  the  buoyancy  equation  with  ^  (t,  z)  as  a  dependent 
variable.  Note  that  (t,z)  is  the  total  upward  mass  flux  in  the  boundary 
layer  integrated  over  the  thickness,  6,  of  the  boundary  layer. 

The  new  form  of  the  buoyancy  equation.  Equation  (60),  can  be  solved  for 

tne  total  mass  flux,  once  the  variables  Acand  aw  are  properly 

determined.  Analoqous  to  the  previous  procedure,  a  and  a„  are  expressed 

w  ^ 

in  terms  of  staircase  functions  in  the  stationary  grid. 


(t,  0  -  (t) 

(t,z)  •  (t) 


for  z,  ,  <  z  <  z, 
k- j.  k 

and  k  -  1 ,2, . . . ,n . 


where  by  definition  '£  .  (t)  and  a  ,  (t)  are 

w,k  C,K 


\,k  =  9  3  (W  -  VO 


\,k  9  3  ^9c,k  -To  ^ 


After  substituting  for  aw  and  ac  in  Equation  (60),  ana  noting  that  they 
are  constant  in  the  interval  z.  ,  <  z  <  z.  ,  the  buoyancy  equation  will  take 

K- 1  K 

the  following  form: 


a  3/4  _  1/3 

Hz  =  CoVk 


for  z(<_1  <  z  <  zk  (61) 


Integration  of  the  modified  buoyancy  equation,  Equation  (61),  results  in 


*3/4(t,z)  =  .1,3/4 


^tlk-l 


(t'  +  Cc/J,'k  (t)(z 


zk-p  for  zk-l 


<  z  <  z. 


and  evaluating  it  at  z  -  for  all  k  =  l,2,...,n  yields 


Equation  (62)  is  the  general  solution  of  the  buoyancy  equation, 

Equation  (61).  It  calculates  the  value  of  the  total  mass  flux,  i.e. 

£>...»  at  any  grid  point  along  the  boundary  layer.  Obviously,  since 
there  is  no  flow  in  the  bottom  of  the  tank,  ^  g  =  0. 

A  more  accurate  form  of  the  bouyancy  equation  may  be  obtained  once  the 

variables  a  and  a.  are  expressed  as  staircase  function  approximations.  The 
w  c 

resulting  form  of  the  bouyancy  equation.  Equation  (61),  would  be  supple¬ 
mented  by  a  j  unp  condition  on  $  at  each  z^,  i.e.,  the  assumption  of  conti¬ 
nuity  of  pj.  at  each  z^  would  be  relaxed.  This  jump  condition  could  be 
derived  as  follows:  Integrate  each  term  of  Equation  (60)  from  z  -  6  to  z, 

K  K 

+  <$  where  6>0.  Then  let  6+0.  The  result  would  be  an  expression  from  which 
the  jump  in  ^  at  can  be  calculated. 

Another  quantity  essential  to  calculation  of  the  core  top  temperature 
is  the  upward  heat  flux  Q^nd(t,  z),  with  Tc(t,  z)  as  its  reference  tempera¬ 
ture,  integrated  over  the  thickness,  6,  of  the  boundary  layer.  It  is 


z) 


C  /  L 
pJo 


(T  -  Tc)  diji 


a  dq 


After  substituting  from  Equation  (59),  the  above  equality  becomes: 

q'bnd(t>z)  =  Cz  \  (t,z)t|/l,(t>z)  (63) 

where 


0.937  vj) /4C 
_ £ _ 2. 

6  g(Pr)l/2 
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At  the  core  top  the  total  mass  and  the  convective  heat  issuing  from  the 
boundary  layer  are  defined  as; 


*t  top  (t,“  *t  <l’  V  <64> 

4bnd  top  CO'  5bnd  (t,  Ht)  (65) 


The  two  quantities  given  by  Equations  (64)  and 

ture  T  ,  (t)  of  the  liquid  flowing  onto  the 

c  top  3 

layer  according  to  the  following  equation 


(65)  determine  the  tempera- 
top  core  from  the  boundary 


T 


c  top 


(t) 


T 

'c 


(t. 


%nd  top 
~^p^t  top 


(t) 

to 


(66) 


The  equations  describing  the  flow  variables  at  the  core  top,  Equations 
(64),  (65),  (66),  written  in  terms  of  their  corresponding  staircase  repre¬ 
sentation,  will  take  the  following  form;  the  total  mass  flux  becomes 


«l». 


t  top 


(t)  = 


(67) 


Similarly,  the  expression  for  the  heat  flux  at  the  top  of  the  boundary 
layer,  written  in  staircase  form  after  substituting  from  Equation  (63)  into 
Equation  (65),  becomes 


(68) 


Following  the  same  procedure,  the  expression  for  the  core  top  temperature  in 
the  staircase  form  becomes 


Tc,  topd  = 


c,n 


(t) 


%nd  top 

^p^t.top 


(t) 

TIT 


(69) 
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F. 


MODELING  OF  VAPORIZATION  RATE 


The  relation  describing  the  vaporization  and  the  mass  transfer  across 
the  gas  liquid  interface  is  given  in  Equation  {22).  This  equation  can  be 
put  in  the  form 


or 


where 


m"  =  a  In  ( 1  + 


b  m 


m"  -  c 


m"  =  a  In  f  1  +  b  + 


b  c 


m"  -  c 


h 


a  = 


=  'g 

p,g  p3g 


<Tq  “  V 


1/3 


g  6  M/3 

k  Ci .  I  — — 2- 

S  lh  1  Vgag 


C  (T  -  T .  ^ 

p.gv  g  w 

*eff 


(70) 


c  = 


ri 

Leff 


£eff  "  L  +  Cp ,£  (Ti  ‘  Tc  top^ 


q"  .A 
Vri  w 


rz 


The  surface  temperature  T ^ ( t )  is  estimated  according  to  t tie  following 
rel ationshi p: 

T.  =  min  [(1/2)(t^  +  Tc  top),  tj. 

The  boiling  temperature  T^  is  determined  by  the  Clausius-Clapeyron  relation 
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and  is  discussed  in  the  next  section. 

Parameters  riot  previously  defined  are 

kg  =  thermal  conductivity  of  the  gas 

C,,  =  dimensionless  constant  for  horizontal  orientation  of  surface 
lh 

8^  =  thermal  expansion  coefficient  of  the  gas 

v  =  kinematic  viscosity  of  the  qas 
g 

a  =  thermal  diffusivity  of  the  gas 
9 

Equation  (70)  is  a  nonlinear  equation  that  needs  to  be  solved  for  iti". 
This  problem  is  a  particular  case  from  a  general  class  of  problems  in  which 
roots  are  to  be  found  for  an  equation  of  the  form 

f(x)  =  0 

Numerical  methods  for  solving  such  problems  are  discussed  extensively  in  the 
literature.  To  convert  Equation  (70)  to  this  form,  transpose  the  right  side 
to  the  left  side  of  this  equation,  and  replace  m"  by  the  symbol  x.  This 
result:  in  the  following  expression  for  the  function  f(x). 

f(x)  =  x  -  a  In  (1  +^~)  (71) 

Then  the  value  x  =  x^,  where  f  =  0  is  the  solution,  m",  of  the 
transcendental  equation,  Equation  (70),  i.e.,  m"~  x^. 

To  develop  a  mjnerical  method  of  solution,  some  properties  of  the 
function  f(x)  need  to  be  determined.  Assume  that  a  >  0,  b  >  0,  c  >  0,  which 
will  normally  be  the  case.  Note  the  properties, 

f(x)  undefined  for  0  <  x  <  c 

f  +  as  x  +  c  from  the  right 

f  -»•  +<»  as  x  >  +® 
df  .  .  , 

•^  >  1  for  c  <  x 
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It  follows  that  there  is  one  solution 


and 


f(x)  <  o  for  c  <  x  <  xQ 
f ( x)  >0  for  xQ  <  x 


An  iterative  scheme*  which  is  a  combination  of  Newton's  method  and  the 
bisection  method  (Reference  29)  is  used  in  the  numerical  solution  of  this 
equation.  For  any  iterate  x,  if  f(x)  >  o,  then  the  next  iterate  is  taken 
half  way  between  the  points  c  and  x,  unless  Newton's  method  gives  a  point 
closer  to  x.  This  procedure  is  continued  until  an  iterate  x  is  reached 
where  f(x)  <  0.  After  reaching  an  iterate  that  gives  f  <  0,  Newton's  method 
is  used  for  calculating  all  the  remaining  iterates.  It  can  be  shown  that  if 
f  <  0  for  one  iterate,  then  all  the  remaining  iterates  from  Newton's  method 
will  also  give  f  <  0  . 


The  input  quantity  f 


the  solution  xQ.  If 
>  1  for  x  >  c. 


tol 
<  f 


is  the  absolute  degree  of  accuracy  desired  in 


tol 


,  then  it  follows  that 


x-x. 


<  f.  .  because 
tol 


G.  MODELING  OF  GAS  AND  THE  UNWETTED  WALL  TEMPERATURE 

The  eouation  for  the  unwet.ted  wall  temperature  T^(t)  was  derived 
earlier,  Equation  (1).  It  satisfies  the  following  ordinary  differential 
equation. 


dT 

p  C  5  — ^  =  n" 

Ui  I 


WWW 


dt 


-  q" 


rro 


q  •  -  q 

M  rri  ^  eg 


(72) 


*  The  method  was  developed  only  to  the  point  where  it  satisfied  the  present 
needs.  Further  investigation  would  be  desirable  to  assess  its  limitations 
and  its  capability.  The  method  of  solving  the  vapori zation  equation  can 
easily  be  improved.  The  bisection  method,  which  was  used  in  part  of  this 
calculation,  is  slow.  The  number  of  iterations  can  be  reduced  by 
modifying  this  part  of  the  calculation. 
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where 


q"  =  ea  T4  -  T4 
Mrro  0  \  wg  « 


q"  =  ea  I  T4  -  T4 
rr  i  \  wg  i 


n  «i  _  r 

qcg  ~ 


(j  „  -  T  V*/3 

\  w9  g^1 


The  coefficient  Cq  was  defined  earlier,  it  reads. 


C  =  k  Cv 
9  9  In 


1/3 


The  gas  temperature  equation,  Equation  (17),  is  modified  to  allow  also 
for  an  open  vent  configuration.  The  gas  temperature  T^  satisfies  the 
following  ordinary  differential  equation: 


dT  ,  N-ir 

— -r^-  =  f  M  C  +  i*i  C  ]  1  q"  A 
ot  ^  a  a,v  v  v,vy  I  Hcg  w 


^A  +  'MCV,P  Ti  *  Cv,v  M 


-  ftQPV/  (Ma  +  Mv) 


(73) 


where 


q"  =  in "l.  ,,  -  q" 
C£  eff  Hr$, 


The  mass  inflow  rate  into  the  ullage  vulume  due  to  the  evaporation 


is 


Mi  =  A£  m" 


(74) 
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The  mass  outflow  rate  M  from  the  ullage  volume  through  the  vent  is 


for  P  <  P 


M  =  /  M  P  -  P  J/aP  for  P  <  P  <  P  +  aP  ,  (75) 

o  \  \  o/max\  vnt/  vnt  vnt  vnt  vnt 


o/max 


for  P  +  aP  ,  s  P 
vnt  vnt 


where 


M  =  AC  .  2  P  -  P  .  M  +  M  1/V 
o/max  d  \  atm/ \  v  a/ 


and  Pvnt  is  the  vent  activation  pressure  such  that  the  vent  is  closed  and 

open  for  P  <  p^  and  P  >  Pvn*.s  respectively.  The  quantity  aP  ,  when  it 

is  set  greater  than  zero,  gives  a  model  of  a  vent  that  is  partly  open  in  tne 

oressure  range  P  .  <  P  <  P  ,+aP  x. 

a  vnt  vnt  vnt 

The  mass  M  (l)  emu  M  (L)  of  vapor  and  air,  respectively,  in  tne  ullage 

V  d 

volume  satisfies  the  following  ordinary  differential  equations 


.  .  IM.-MM  /  K  +  M 
d  t  /  1  o  v  \  v  a 


for  P  <  P 


for  P  ,<  P 
vnt 


'-w( 


for  P  <  P 


for  Pvnt<  P 


The  four  simultaneous  ordinary  differential  equations,  numbered  (72),  (73), 
(76),  and  (77),  are  solved  by  a  fourth-order  Runge  Kutta  numerical  method 
( Reference  27) . 
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H.  MODELING  OF  PRESSURE  AND  BOILING  TEMPERATURE 


According  to  Equations  (18)  and  (20),  the  total  pressure,  P(t),  of  the 
gas  is 


P(t) 


Ma(t)  Mv(t)  ]RTg(t) 
TmoTJ^  +  TmoTJ^J  V~ 


(78) 


The  boiling  temperature  T^(t)  of  the  liquid  at  the  total  pressure  P(t) 
is  given  by  the  Clausius-Clapeyron  equation. 


r 


Tb(t)  - 


T, 


(mol)vLv 


in 


?(t) 


(79) 


where  Px  is  the  vapor  pressure  at  the  reference  temperature  Tt;  i.e.,  Tx  is 
the  boiling  temperature  when  the  total  pressure  is  P^ 


I.  INITIAL  VALUE  CALCULATION 


The  initial  values  of  the  physical  parameters  in  the  ullage  volume. 

i.e.,  vapor  pressure  P  (0),  vapor  mass  M  (0),  and  the  amount  of  air  M  (0) 

v  v  a 

are  calculated  using  Cl ausius-Cl  apeyron  and  ideal  gas  law  relations. 


Using  the  Cl  ausi us-Cl apeyron  equation,  the  initial  vapor  pressure  in 
th-:  tank  is  obtained  from  the  following  relation; 


Pv(0) 


Pt  exp 


(mol)  L 

V  V 

R 


(80) 


The  corresponding  initial  mass  My(0)  of  the  vapor  in  the  tank  is  calculated 
by  substituting  Py(0)  from  Equation  (80)  in  the  ideal  gas  law  relation; 


(mol)  V  P  (0) 

M  (0)  = - i - l - 

V  Rig  (0) 


(81) 


i 

I 
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Assuming  that  the  initial  total  pressure  inside  the  tank  is  equal  to  the 

ambient  atmospheric  pressure  P  ,  the  initial  mass  MjO)  of  the  air  is 

a  cm  a 

obtained  front: 


M  AO) 

a 


-  pv(o). 
RTg  (0) 


(82) 


J.  DESCRIPTION  OF  THE  CODE 


A  modular  code  using  the  FORTRAN  77  language  was  written  to  numerically 
solve  the  governing  equations.  The  code  was  developed  in  the  VUE  (VOS  UNIX'" 
Environment)  operating  system  on  a  Harris  800  computer.  The  code  consists 
of  one  main  program  and  seventeen  subroutines.  A  description  of  each  of 
these  program  units  is  presented  below.  A  copy  of  the  source  code  is  given 
i  n  Appendix  B. 

Program  TANK  is  the  main  program.  To  allow  flexibility,  tne  program 
has  two  time  variables,  and  t^,  for  gas  and  liquid,  respectively.  Tne 
two  correspond  ing  time  increments  can  be  related  so  that  At^  =  m  At^,  where 
m  is  a  positive  integer.  However,  in  this  calculation  m  =  1  was  used. 

The  main  program  performs  the  following: 

1.  Initially,  it  calls  the  INPUT1,  INPUT?,  INPUT3,  and  INITVAL 
subrouti nes . 

2.  It  carries  out  the  iterative  nunerical  procedure  by  calling  GKUTTA 
for  each  updating  of  the  gas  variables  by  a  time  increment,  ,  and  by 
calling  the  BUOY,  CORE,  LKUTTA,  and  RE2C0RE  subroutines  for  each  updating  of 
the  liquid  variables  by  a  time  increment.  At^. 

3.  It  controls  the  output  time  by  the  variable  t  .  When  t  becomes 

pr  i 

1 arger  than  t  ,  the  program  calls  WL  and  DVG  to  complete  the  calculations. 
It  then  outputs  the  desired  liquid  variables  as  well  as  the  desired  gas 
variables  on  the  VOS  files  called  OUTL  and  OUTG,  respecti vely.  Similarly, 
it  outputs  the  interface  and  the  remaining  variables  on  the  VOS  file  OUT  I . 
Finally,  it  increments  tpr  by  the  amount  Atpr  to  determine  the  next  output 
time. 
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4.  The  program  stops  when  t  exceeds  t 

pr  max 

Subroutine  INPUT1  assigns  and  calculates  constants  and  initial  values. 

The  initial  values  assigned  in  this  subroutine  are  the  initial  time,  t  = 

o 

TIMZ,  and  the  initial  temperature,  T  =  TZ.  To  change  the  values  to  be 
assigned  in  INPUT1,  the  source  file  which  has  this  subroutine  is  edited  and 
recomp i 1 ed . 

Subroutine  INPUT?  reads  the  VOS  file  INTK  to  get  the  following  values: 

the  maximum  time  t  -  TIMMX;  the  time  increments  At  =  DTIMl,  At  =  DTIMG, 
max  e  g 

and  At.  =  OTIMPR;  the  parameters  T.,.  -  TDFTOL,  x,  ,  =  XTOL, 

pr  r  dftoa  tol 

f  j  =  FTOL,  and  n  ^  =  NTOL  which  determine  the  accuracy  of  the  iterative 

calculation  in  Subroutine  EVNEWT;  the  number  n  =  N  of  computational  cells  in 

the  stationary  qrid  z  ,  z,,...,z  ;  and  the  vent  parameters  V  ,  =  VOPT  = 

o  i  n  opt 

VENTOPT  =  ("Y"  for  vented  option,  "N"  for  unvented  option),  P  =  PVNT,  and 
APvnt  =  DPVNT,  It  then  outputs  all  these  input  values  to  the  V OS  file 
OUTTK. 

Subroutine  INPLIT3  interact i vely  obtains  either  a  fire  diameter,  DIAMF, 
which  it  converts  to  heat  flux  by  a  table  (when  provided),  or  the  heat  flux, 
q'jt  =  QF,  directly.  It  then  outputs  these  input  values  to  the  terminal  an¬ 
al  so  to  the  VOS  file  OUTTK. 

Subroutine  INITVAL  initializes  the  temperature  variables  T.  ,  T 

wb*  ws,k* 

T  ,,  T  ,  T  and  T  ,  .  It  also  intializes  the  ts\o  time  variables  t  and  t  . 
c,k  wg  g  ct  19 

It  calls  INI TM  to  calculate  initial  values  for  some  gas  variables.  It 

calculates  the  stationary  vertical  grid  points  z0,  z1,...,zr)  for  the  liquid 
part  of  the  system. 

Subroutine  INI TM  calculates  the  initial  values  for  the  vapor  pressure 
Pv ,  for  the  mass  of  vapor  in  the  ullage  volume,  and  for  the  mass  of 
air  in  the  ullage  volume. 

Subroutine  GKUTTA  updates  the  gas  variables  T  ,  T  ,  M  and  M  over  a 

wg’  g*  v  a 

time  increment  At  by  using  a  fourth-order  Runge  Kutta  method  for  solving 
the  system  of  simultaneous  ordinary  differential  equations  for  the  gas 
variables.  It  obtains  the  first  derivative  with  respect  to  time  of  the  gas 
variables  by  calling  subroutine  OVG. 
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Subroutine  OVG  calculates  PK ,  Mq,  and  the  first  derivative  with  respect 

to  time  of  the  gas  variables  T  ,  T  ,  M  and  M  .  It  calls  subroutine  WG  to 

wg’  g  v  a 

evaluate  other  variables  appearing  in  the  differential  equations  of  the  gas 

variables.  The  other  variables  are  functionally  dependent  on  the  gas 

vari ables  and  on  T  .  . 

ct 

Subroutine  WG  calculates  the  variables  which  are  functionally  dependent 

on  the  qas  variables  T  ,  T  ,  M  ,  M  and  also  on  T  . .  First,  it  calculates 
3  wg  *  g*  v’  a  ct  * 

q"  ,  q"  ,  F*  P  ,  P,  and  T,  .  It  calls  subroutine  SURFTEMP  to  obtain  the 
Vro  eg  a’  v’  ’  b  -  ■ 

interface  surface  temperature  .  It  then  calculates  q^ri  ,  q^,  ana  Le^f. 

It  calls  subroutine  EVAP  to  obtain  the  value  of  the  evaporation  rate  m". 

Then  it  calculates  q 

Subroutine  EVAP  calculates  the  value  of  the  evaporation  rate  in".  It 
calls  subroutine  EVNEWT  to  solve  the  transcendental  equation  for  m" . 

Subroutine  EVNEWT  numer ical ly  solves  the  transcendental  equation  for 
the  evaporation  rate  m"  using  a  combination  of  Newton's  method  and  the 
bisection  method. 

Subroutine  SURFTEMP  gives  an  estimate  of  the  value  of  the  surface 
temperature  T  . 

Subroutine  BUOY  calculates  the  boundary  layer  mass  flux 

:b  ,  (t  )  =  Or  (t  ,  2,  )  for  k  =  0,1,..., n  at  the  stationary  grid  heights  z  , 
t.,k  i,  t  i  k  o 

Zj,...,z  by  using  the  analytical  solution  of  the  discretized  buoyancy 

equation  in  which  the  side  wall  temperature  T  (t,  z)  and  the  liquid  core 

temperature  Tc(t,  z)  are  expressed  as  staircase  functions  in  the  stationary 

grid.  It  then  assigns  the  value  for  the  boundary  layer  mass  flux  at  the  top 

iii.  .  (t  )  =  iii.  (t  ).  It  calculates  the  value  for  the  boundary  layer 

yt  top  v  vt,n  v 

heat  flux  at  the  top, 


^bnd  top^j^  ^bnd  i* 


Subroutine  CORE  calculates  tne  grid  heights  z^,  z[  ...,z^  at  time 

t  +  At^  of  the  moving  grid  which  initially  coincides  with  the  stationary 

grid  at  time  t  and  moves  with  the  liquid  core  during  the  time  step  At  .  It 
£  t 

then  sets  z'  ,  -  z  .  The  tenperature  T  .  (t  +  At  )  at  the  liauid  core 
n+1  n  r  c  top  li 

top  is  calculated.  The  temperature  T  (t  +  At  ,  z)  of  the  liquid  core  at 

C  %  2, 

the  end  of  the  time  step  At.  is  calculated  as  a  staircase  function  in  the 

i 

moving  grid,  i.e. 

T  (t  +  At  ,  z)  =  T'  .  (t  +  At.  )  for  z,1  ,  <  z  <  z' 
cv  i  t  c,k  v  a  i'  k-1  k 

for  k  =  1,2, . ..  ,n+.l. 

Subroutine  LKUTTA  updates  the  variables  T  ,  and  T  .  for  k  =  l,...,n 

wb  ws,k 

by  a  time  increment  At  using  a  fourth-order  Runge  Kutta  method  to  solve  the 
system  of  ordinary  differential  equations  for  the  variables.  It  obtains  the 
first  derivative  with  respect  to  time  of  these  variables  by  calling  subrou¬ 
tine  DVL. 

Subroutine  OVL  calculates  the  values  of  the  first  derivative  with 

respect  to  time  of  the  wetted  wall  temperatures  T  ,  and  T  ,  for  k  = 

wb  ws,k 

l,2,...,n.  It  calls  subroutine  Wl  to  evaluate  other  variables  appearing  in 
the  differential  equations  for  the  wetted  wall  temperatures.  The  other 
variables  are  functionally  dependent  on  the  wetted  wall  tenperatures  and 
also  on  the  core  temperature  T  ,  for  k  =  l,2,...,n. 

C  9  K 

Subroutine  WL  calculates  the  variables  q^rob9  q^,  ^rros  k’  anci 
q"  v  which  are  functionally  dependent  on  the  wetted  wall  temperatures 

C  jL  o  i  K  j 

^wb’  Tws  k  and  core  temperature  Tc  ^or  k  =  l,2,...,n. 

Subroutine  REZCORE  takes  the  following  staircase  representation  of  the 
core  temperature  in  a  moving  grid 

T  (t  +  At  ,  z)  =  T'  ,  (t  +  At  )  for  z\  .  <  z  <  z/ 

C  £  t  C,K  £  £  k-1  k 

and  converts  it  into  the  following  staircase  function  in  the  stationary 
grid 


Tc(t«  + 


it 


z’ 


z)  '  Tc,k  +  “V 


for  zt_j  < 


z  <  z, 


56 


SECTION  V 
NUMERICAL  RESULTS 


The  physical  processes  occurring  inside  a  fire-exposed  tank  containing 
liquid  fl  ammabl  es  were  identified  and  analyzed,  and  the  corresponding 
governing  thermohydrodynamic  equations  describing  the  physical  phenomena 
were  derived  in  Section  III.  Numerical  solutions  of  the  governing  equations 
and  modeling  of  the  tank  response  prediction  to  an  external  thermal  lead 
were  presented  in  Section  IV.  This  section  describes  a  special  numerical 
scheme  that  was  developed  to  calculate  the  response  of  the  individual  com¬ 
ponents  comprising  the  fire  tank  system.  An  integrated  response  methodology 
was  developed  based  on  the  response  of  wetted  and  unwetted  wall  tempera¬ 
tures,  gas  and  liquid  phases,  convective  turbulent  boundary  layer,  interface 
mass  evaporation,  and  the  overall  tank  pressure  rise. 

Calculations  are  made  to  predict  the  response  of  a  vented  and  an 
invented  55-gallon  (206-liter),  18-gage  steel  drum  to  a  uniform  external 
heating  of  7  w/ cm2 .  For  the  vented  cases,  the  vent  diameter  is  7.6  mm  and 
the  vent  opening  pressure  is  41  kPa  gage.  The  vent  coefficient  or  discharge 
used  in  the  calculations  is  1.0  for  the  Freon™  tanks  and  0.25  for  the  water 
tanks.  The  drums  are  filled  up  to  95  percent  filling  capacity  either  with 
water  or  Freon-113'".  In  the  cal  cul  ati  ons ,  it  is  assumed  that  the  drums  are 
completely  engulfed  and  the  total  incident  heat  consists  of  radiation  as 
well  as  convection. 

The  response  predictions  of  the  wetted  side  wall  and  bottom  wall  tem¬ 
peratures  for  vented  and  (invented  configurations  are  shown  in  Figures  9  and 
10  and  11  and  12,  respectively.  The  calculation  shows  that  the  temperature 
of  the  wetted  walls  strongly  depends  on  the  liquid  in  contact  and  is  not  at 
ail  affected  by  whether  the  tank  is  vented  or  unvented.  It  is  also  noted 
that  the  bottom  wall  invariably  shows  higher  temperature  than  the  side 
wa'ls.  Unlike  the  wetted  walls,  as  might  he  expected,  the  temperature  of 
the  dry  walls  in  the  ullage  volume  is  higher  than  the  wetted  wall 
temperature  and  remains  the  same  for  both  liquids  (Figures  13  and  14). 

The  predicted  mass  and  heat  fluxes  for  water  and  Freon-113™  in  the 
boundary  layer  for  vented  and  unvented  lank  configurations  are  shown  in 
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igure  9.  Wetted  Sidewall  Temperature  Time  Profile  for  Closed  Vent  Tank  Configuration 


Figure  10.  Wetted  Sidewall  Temperature  Time  Profile  for  Vented  Tank  Configuration 
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Bottom  Wall  Temperature  Time  Profile  for  Closed  Tank  Configuration 


Figure  12.  Bottom  Wall  Temperature  Time  Profile  for  Vented  Tank  Configuration 


Figure  13.  Dry  Wall  Temperature  Time  Profile  for  Closed  Vent  Configuration 
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Figure  15.  Boundary  Layer  Hr at  Flux  Time  Profile  for  Closed  Vent  Configuration 
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Figure  16.  Boundary  Layer  Heat  Flux  Time  Profile  for  Vented  Configuration 
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guration 


Figures  15  and  16  and  17  and  18,  respectively.  Figures  15  and  17  show  that 
the  rate  at  which  heat  and  mass  are  transported  along  the  boundary  layer  is 
independent  of  the  tank  configuration  and  is  controlled  by  the  physical 
properties  of  the  lading.  It  is  noted  that  the  Freon-113”  high  rate  of  mass 
flux  in  comparison  with  water  is  evidence  of  this  phenomenon. 

The  predicted  pressure  time  profiles  in  the  ullage  volume  are  shown  in 
Figures  19  and  20.  As  shown  in  Figure  19,  in  an  unvented  tank,  the  inside 
pressure  rises  quickly  (in  less  than  1  minute)  to  a  level  that  no  ordinary 
container  can  withstand,  and  consequently  the  tank  ruptures.  However,  the 
pressure  rise  inside  a  vented  tank  is  more  manageable.  The  pressure  rises 
sharply  until  it  reaches  the  vent  operating  pressure  value  at  which  venting 
takes  place  (Figure  20).  The  oscillation  observed  is  caused  by  the  opening 
and  closing  of  the  vent  before  it  becomes  stable.  Beyond  this  point  the 
pressure  continues  to  rise,  hut  at  a  slower  rate  because  of  the  expulsion  of 
the  tank  contents.  With  increasing  time,  the  internal  pressure  rises 
steadily,  causing  the  tank  to  rupture  in  less  than  200  seconds.  The  dry 
wal'  temperature  for  hoth  Freon”  and  water  was  about  900  K  for  the  closed 
vent  tanks  (Figure  13).  For  the  vented  tanks,  the  dry  wall  temperature  was 
essentially  the  same  as  for  the  closed  vent  tanks  for  Freon-113”  and 
slightly  less  for  water  (Figure  14).  The  wetted  sidewall  temperatures  were 
essentially  the  same  for  both  the  closed  vent  and  vented  tanks;  that  is, 
about  400  K  for  water  and  500  K  for  Freon”  (Figures  9  and  10).  The  model 
prediction  of  other  physical  parameters  sucn  as  mass  of  vapor  in  the  ullage 
volume,  hoi  ling  temperature,  and  gas  temperature  are  shown  in  Figures  21- 
26. 


The  predictions  from  the  calculations  shown  in  Figures  9-25  are  off  as 
a  result  of  inadvertent  omission  of  a  constant.  The  liquid  density  factor 
(RH0Z)  was  omitted  from  the  expression  for  the  constant  CZ  in  both  the  water 
and  Freon”  versions  of  subroutine  INPUT1.  This  resulted  in  very  small  pre¬ 
dicted  boundary  layer-  mass  and  heat  fluxes.  The  code  was  corrected  and  was 
also  improved  to  include  the  effect  of  venting  on  the  rate  of  change  of  the 
gas  temperate ~e.  The  corrected  code  predicted  a  boundary  layer  mass  flux 
that  appears  to  be  much  too  large,  because  of  project  resource  limitations, 
this  was  not  investigated  further. 
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figure  21.  Mass  of  Vapor  in  Ullage  Volume  Time  Profile  for  Closed  Vent  Configuration 


i gu re  22.  Mass  of  Vapor  in  Ullage  Volume  Time  Profile  for  Vented  Configuration 


Figure  23.  Boi 


700 


Figure  26.  Gas  Temperature  Time  Profile  for  Vented  Configuration 


SECTION  VI 
CONCLUSIONS 


The  governing  equations  describing  the  interactive  processes  occurring 
between  the  fire  and  the  tank  have  been  solved  numerically,  and  1  model 
predicting  the  thermo* luid  procesres  leading  to  occurrence  of  8L.EVE  ha3  been 
developed.  The  calculations  were  made  for  vented  and  unvented  tanks 
containing  either  water  or  Freon-113'".  The  prediction  demonstrated  internal 
consistency  amongst  the  parameters  indicative  of  the  moaei  performance . 

Based  or-  ?  1  y s i s  of  the  calculate1  prediction,  the  following  general 

observations  made- 

1.  The  pressure  rise  in  tire  ullage  volume  and  occurrence  of  BlEVt 
strong  y  depend  on  the  tana's  initial  filling  density. 

2.  The  physical  properties  of  lading,  especially  the  heat  of 
vaporization,  s  ;  gni ficantly  contribute  to  the  internal  pressure  rise. 

3.  The  liquid/gas  interface  is  the  most  thermally  active  region,  and 
a  decreasing  temperature  gradient  exists  in  tne  downward  direction. 

<■'. .  Only  the  liquid  boiling  temperature  is  affected  by  configuration 
(vented  nr  unvented)  of  the  tank;  all  other  temperatures  remain  unaffected, 
regardles  of  the  tank  configuration. 

Specif  'ally,  to  delineate  several  fundamental  issues  concerning  this 

development,  'hi  following  preuiction  case  is  considered.  The  calculated 

response  for  a  vented  55-gallon,  18-gage  steel  drum,  engulfed  in  an  intense 

fire  (7.0  )  and  filled  to  95  percent  of  its  capacity  with  water  showed 

cm2 

that 

1,  The  tank  failure  pressure  is  reached  at  about  £00  seconds.  The 
calculated  pressure  for  a  closed  tank  approached  values  much  higher  than  for 
the  vented  tank  for  the  same  time. 


1.  The  dry  wall  and  wetted  wall  temperatures  are  changed  very  little 
by  tank  venting. 

3.  The  internal  pressure  in  a  vented  tank  rises  much  more  slowly  tha^i 
the  pressure  in  a  closed  tank. 

Similar  calculations  were  made  using  Freon-113'"  in  place  of  water  while 
maintaining  the  other  variables.  As  was  expected,  the  calculated  response 
predictions  for  Freon-1131"  were  quite  different  when  compared  with  tne  case 
for  water. 

The  corrected  code  predicts  a  significantly  higher  boundary  layer  mass 
flux.  Causes  of  this  problem  might  be  the  simplifying  assumptions  made  and 
the  empirical  fits  used  in  the  derivation  of  t he  model  equations,  or  the 
numerical  scheme  used  to  solve  the  equations.  Project  resource  limitations 
'id  not  permit  determination  of  the  cause  and  the  resulting  effects  on  the 
predictions. 

The  failure  response  of  the  tank,  and  ultimately,  the  occurence  of 
I3LFVF  also  depend  on  the  size,  and  the  mechanical,  thermal  and  physical 
integrity  of  the  tank.  Those,  together  with  the  variations  in  heat  fluxes 
from  various  fire  sizes,  as  well  as  differences  in  the  heat  input  caused  bj. 
the  relative  location  of  the  tank  with  respect  to  the  fire,  indicate  that 
every  response  prediction  calculation  requires  a  myriad  of  different  inputs. 
This  becomes  a  laborious  and  time-consuming  operation,  not  consistent  with 
the  purpose  for  which  the  model  was  developed.  The  mission  of  this 
development  is  to  quickly  provide  to  the  firefighters  the  essential 
predictions  and  the  risk  assessment  data  necessary  for  safe  operations  and 
successful  extinguishment  of  a  fire. 

To  achieve  this  goal  without  compromising  accuracy  of  prediction,  it  is 
recommended,  once  the  model  is  refined  and  verified,  that  prediction 
calculations  be  made  for  most  probable  ranges  of  fire  tank  situations  and 
that,  a  Bl.FVE  prediction  data  base  be  created.  The  results  from  this  data 
base  could  then  be  plotted  into  a  series  of  reference  charts,  or  stored  in  a 
computer  and  selected  arid  displayed  on  the  screen. 
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APPENDIX  A 


DERIVATION  OF  EQUATION  OF  MASS  TRANSFER 
STAGNATION  FILM  THEORY 


Consider  a  th'n  stagnant  layer  of  fluid  with  thickness  6,  Figure  A-l 
(Reference  14,  p .  135) . 


1 


AY 


-*  X 


Figure  A-l.  Stagnant  Film  Layer 

Assune  that  due  to  the  high  temperature  gradient,  mass  is  vaporizing  at 
the  free  surface.  Define  the  layer  thickness  so  that 


K/6  *  h 


(A-l) 


The  conduction  equation  in  this  layer  is 


2 

r  ..  3T  _  .  3  T 

‘V  w  -K  — 

3^ 


( A-2 ) 


If  m"  is  the  mass  rate  of  vaporization  per  unit  surface  area,  the  conduction 
equation  (A-2)  can  be  modified  using,  v  =  to  read 


Cm"  il  =  K  -U 
p  sy 

The  differential  equation,  Equation  (A-3),  is  integrated.  We  have 


(A- 3) 


( A-4 ) 
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Applying  the  following  boundary  condition,  one  obtains 


at  y  *  o,  T  *  Ti  ;  at  y  =  6,  T  =  TQ, 


The  solution  then  becomes 


T a  '  T i  =  C1  —  E*P  (Cnm"  6/k)  -1  ( A— 5 ) 

y  1  C  ft"  p 

p 


The  energy  balance  at  the  surface  gives 


m"  Leff  =  q"con 


6"  a  k  -  +  n" 

f'.A  a  yfo  ■  r,i 


From  equation  (A-4) 


„  =  C,K 
y=o  1 


and  hence  frcm  Equation  ( A- 6 )  we  have 


(A-6) 


CjK  in  Loff  -  q 


An 


eff  (  r,z 


(A-7 


Substituting  for  C,K  from  Equation  (A-7)  into  Equation  (A-5)  one  obtains 


CP  (T9  '  T0 


Leff  '  q  r,ft 
m"  / 


exp  (  Cm"  6/k  (  -1 


(A-8) 


Equation  (A-8)  describes  the  mass  transfer  rate  which  can  be  solved  for  ;ii" 
to  give  Equation  (22) . 


APPENDIX  B 


PROGRAM  LISTING 
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.  <Ln.  C  (July  -  •„  •  1  -  i  ) 

T  fh  ‘olio  n  in;  t  -•  )  -  r  f  hi  '  r.  1  v  o  i  lha  d  £  r  1 1  ?  1  Jltf-r»n1jrl 

l;  U:r  ;  ci  stf  i5ir'  the  tt  trschycricyn?  rit  ur:;>;;’5  e'eurrme 
in.u;  *  t‘n>  r:-;i,iUnv  rrcii  t  l  r  ?  itpiorantp'.,  Tna  n  r  c  5  r  ?  m  is 
liquids  5rcu  jii  t :  cllou  >;tr?r  njet  t  i  l  x  cr  fire  cis»et»p  to  b  a 
x  n  p  u  x  • 

Ini  :  :n  -  r.fr  ;  ir;l;,r  in  thi-  r OPT  PAN  7?  1  ;  n  ■;  l c  e  .  T  h  1  code 
•  i:  oavtict;."1  in  xh?  v  „  i  (VCS  UNIX  Environ n. ant)  ap«r?tinc  system 
cn  c  c:ifut>r.  It  ccn-ist  s  ct  one  -rein  p»-ocrpii  and 

sjvsntiin  iuorou tires. 

«»*»*»  .»«••>««  ...... t**.********.********************** 


i  i.  :  '  C  l  T  * i  *IiPlT  1 
C  rlKii  Version 

CC! VON  /rnYSCM/  5/G/c£TP 
C  ON  /GcO'CM/  -L/ XLL/ I-L  /-» /Ac  /  v,  VT 

C  C  I-  ,v  C  M  /  ;•: i  1  ?  P P CN /  *!'Clii/>.MClV,Vne/flnOZx?cT£/XLV,XLxT1,P1 
C  Cfy-Z-ti  /-clTCPCN/  CP/CP  L/C  VP/CW.Ci>VxCPG 
CCKM0.N  /.-C£TPlC«/  KCu-xEPSlxCZxCc/C£/CS/CC-/Ci  V 
C  L  .v  v  0  N  /  I  M  T  C  M/  T > Z, TZ /„ ?, T IN? 

L  ON  /  V  c f  T  C  •'*  /  MM/CjVM/PVM,OFVM 
«  1  i  .  3  C  -  £  7 

j  -  7.1 

fIIi'1  =  1  .  L 1  p  2  ;  i 

fi  L  i  0 . 7  1  -  1 

-  L  =  0 

XLL  -  k  :T(  J  .  lil  ;  ;»£L) 

-  *  -  C  .  <<  3  j  7 

-  :  -  £  l 

»•  *  it-. 

V  T  -  £  l  *  f-  L 
X  e  L  £  z  £  >  .  u 
( ■’■  C  L  v  *  1  i  •  C 
=  1  .  i  -  7 
*  -  -  1  ;  n . : 

a  i_  ~  1  i 

A  L  V  *  A  . 

T 1  *  it..; 

-  1  =  1  .7 1  -  -  3 

l  (■  .  =  l  . r  •:  i 

cp  -  Cfl 

v.  1i  f  -  1.L*.6oxj 

C  VV  =■  1  .  C  L  e  :  3 
Lta  =  7  ,  ?  5  ;  c  2 

C  Pu  -  C  V  r 
i  C  C  n  -  4  .  :  < ; J 
i  P  =  0 .  V 
i 1  -  5  .  c  ?r“ t 
cfil  -  c  ?  «  S  I 
XNLC  =■  4.  .3- 7 
pl'  s  o .  3 

C  i  -  °rCZ*C?*C  x*iuC*atTi«:-/((:.?**4.)*PK**<)  )**(1  .C/3.C) 
x  ixC  -  7 . 5 3  £  *  * 
i  ;  T  AG  =  i.li-i 
Xf.L'C-  -  1.4? 

<\  ri  C  t  L  P  h  -  7.7 
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CALPH  =  0. 71 1 c3 
A  LPHG  1  X  K  C  / (RhCALFH*CAL°H) 

ClhT  =  C.C6 

C 1 ^ I  =  C  .0  = 

Cb  =  XKG*C11^T-C  3=7 AG*G/ ( XNGG*AL PHG )  >**C1. 0/3.0 

CrV  =  XxC-*C1HI*(  C£cT AG*G/ ( XNUC* ALFHG) ) **0 .0/3 .0)  )/CPG 

XNLM  =  C.93?*(5.3)**(4.C/S.G) 

C  f  A  C 1  =  (XMJC**(1  .C/1  2.C) ) / <XNUM*PR**(1 .0/6.0) ) 

CFAC2  -  0. 937*(XNUC**  (1  .0/ A  .0) )  *.CP/ (PR**  <0. 5)  > 

c:  =  cfaci*phc: 

Cl  -  CrACi/ (3=TA*G> 

CIVS  =  C.11 
C  1  P  3  =  C.C* 

C 1  =  CS*C1-=/C1VS 

T  IV  G  =  C  .  0 

T  I  =  2  5  3  .  C 

if  -  7.0  =  4 

Tlx?  =  2  >  5  .  C 

A  V  KT  =  4  ,  5 e 1  =  -5 

COVNT  =  I .0 

PVNT  =  1  . 4  2  fc9; 5 

C  P V  NT  =  G  .  C 

RETURN 

-ML 


SG2RCU7IH:  I M  P  U 1  1 
C  kll:5  Vjrsicn 

C  CI'MCN  /PHY  SC  M/  R/C-/PATM 

COMMON  /G5CMCM/  HL/XLL/ AL/AV./A3/V,VT 

CCFVCN  /MATrPPCM/  X MOL  A / XMC L V, FhC/ 4  “  G  Z  / t = TA / XLV, XLxTI , PI 

CCKMCN  /HSATCPCf'/  CP^CPL/CVP/CVV/CAVxCPC- 

COMMON  /rEATPLCH/  fiC0Ks£PSI/CZ/C2sC$sCS/CGsCEV 

COMMON  /IMTCM/  TIMZ/TZ /tryTINr 

COMMON  /VcNTCK/  AVM/COVNT/PVM/OPVM 

«  =  5.2C9E3 

G  =  9  .  t 

PAT, vi  =  1.0132  =  3 
n  L  "  C . 7191 
a L  =  0.2  333  = 

XU.  =  Z  .C*3Cct(3.14159*AL) 

av.  -  0 •  4  3 2 £ 

it  =  A  l 

V  =  C.CtPiii 

VT  =  AL*^L 

XMC  LA  —  2  9 . C 

XMC  L  V  -  1S.C 

RrlCZ  =  1.033 

kHC  =  Rt-OZ 

=  ET A  =  1  .23-4 

XL  =  2.2533= 

XLV  =  XL 
T 1  =  37  3.0 
PI  =  1.C132E5 
GPL  =  4,1353 
C?  =  CFL 
CVf  =  1  . 9  OE  5 
CVV  =  1.445* 

CAV  =  7.952=2 
CPb  =  CVP 
R  C  C  K  =  4.32=3 
i F  -  0.9 
SI  =  5.475-2 
cPSI  =  E?*SI 
XNLC  -  9. 3c” 7 
PR  =  6.2 

CS  =  RMCZ*CP*(  XNUC*e  5TA*G/ (  <  5 . 3**4  )  «  PR*  *2  )  )**(1.C/2.C) 
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XXG  =  C.cOI 
j  c  T  t  G  -  1 » 1 1 
XfsLC-  = 

f.nCALPt-  -  1.36 
C  -L  Pr-  =  CV? 

F  L  F  rt  G  = 

C  1  h'T  =  C.G  i 
C  1 X  =  C.Ge 

CG  =  XKC*C1nT*(  5ETiG*G/ CXMJG*ALPHu)  )**(1.C/3.C) 

Civ  =  X*6«C1HI»(  (ceT6G*G/(XNLC’'i,wPHCO)**(1  .0/1.0)  )/CPG 
XM..V  =  C.^37«(5.5)»«(i.C/3.C) 

C  f-  £  C  1  -  (XKGC**  (1  .  C  /  1  2.0  )/,(XNsjr-P(.**(1  .0/6.0)  ) 

CF5C2  =  C  .537*  (XNLC**  (1  .0/*.  .C)  )*CP/ (PR**  (0. 5) ) 

CZ  =  CFtC1*RH-JI 

C  2  =  CPSCc/<S5TA*-;> 

civs  =  C.11 

ClrS  =  C .03 

C  s  =  C$*CV-  =  /C1 VS 

T I r  2  =  C.C 

1C  =  25  3.  v 

i.  ■“  -  7  ,  C  ;  4 

TIN.-  =  1 5  l .  v, 

-  vM  =  <..5*15-5 
CCVnT  =  O.cj 
V  VK  7  =  1 5 

0  r  V  -\  T  -  .  C 

\  1  L  v  3 

i  N  C 


i  0  i  ?  C  'J 1  i  X  i  I  '•  P  u  T  2 

CC.V  '"ON  /  T  I  v C  'i /  T  I, vs  X/  Cl  CMC/  07  ZrC  /  Cl  :*.?P/T  ■  , TI«C/  T  I  MPP 

C-jyy  C'i  /  T  C  L  C  •'•'  /  ToFTCU/XTOl/FTOL/NTCl 
COi'^C-'*  / C  v  /  N 

CCv"C*-‘  /V:NTC»/  4V(.  I/CuVNT/PVNT/OPVM 
CiuPaCT:-.  *1  ViNTCrT 
Cr.*-RiCli?  *3  VCTT 


C  r  2  f .  (  5 1  / 

f  IL3  =  '  2GCCCr!Af.  *IMK* ) 

RXCOlii/ 

»  ) 

^s:(55/ 

*  )  T  :  y  (/  X 

■>c20(53/ 

*) 

k-cC<35/ 

*) 

*5  OTIM./  C 

Tirt/  c t : m r s 

■  z  r r  /  ;  c 

*  u  *»  -  v  >  w  / 

*) 

R  cPO ( 5:  / 

«  5 

R  t  C  (  3  5  / 

*1  TOFTCL/ 

XTCl/  F  T  C  L  /  KTCl 

R  cZC(55/ 

*) 

kcXCCSS/ 

*  ) 

S  i  A  C  C  5  5  / 

*  )  5i 

S ;fiC ( 55/ 

»  ) 

Rt::(55/ 

*) 

R  C  A  :  (  5  5  / 

*)  V iNT  OPT  / 

FVNT  /  OFV.'.T 

(  VsNTCPT  .  30.  '  "N"  ) 

GO  TO  1C 

IrCV;NTC?T  .  ;C.  "n") 

GC  TO  1C 

VCFT  =  " 
0  0  T  -J  15 

Hi" 

p  V  N  T  5  1 

.Czii 

V  C  r  T  =  ” 
CCMiM: 

NC" 

CPcN(65/ 

EIL?='2CCCC 

fiR'O'JTTx  '1 

wRITi  Ui 

/  *)  " 

T ;  v  *•  x" 

RPXTcUC 

/  *(:X/c1}.5 

)  '  )  T  I  (■'  K  X 

«  R  :  T  c  ( t  5 

/  *) 

a  R  i  T  z  C  e  5 

/  *)  " 

0  T  !  v  L 

WRITeCci 

/  '  ( ! C  3  x, Cl  3 

.51)')  0  7  I N  L /  cr 

■o  R  i  T  L  ( t : 

/  *) 

Vi  R  :  7  5  (  t  i 

/  '(  " 

TCFTCl 

GTIFG 
C  TING/'  DT1.VP? 

XTCL 


T IFOR" 

F  TCI 
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i  NT  C  L"  )  * ) 

wPiTcltE,  '  C  (  ?X,  1 1 ; . i } / iX,  I:  )  ')  TCFTLL/  XTCL/  F  TGL  /  NT  CL 

»  K  IT  t  (  t  l  ,  *  ) 

hiil  TcUi,  »)  "  N " 

WRITE  (c  E  /  *  (3x,IS  )  ' ) 

hnlTct';/  *  ) 

*)  "  V  EMC  FT  P  V  NT  r  PVNT " 

WRITE (  c  S  r  '  (1  ‘X/ 4  i/2  (  *X, ~ 1 } . S  )  )  '  )  VGPT/PVNT/DPVNT 
WAITi(e5/  *) 

n  C  T  L  K  N 

ENC 

S  0  c  R  u  U  T  IN:  I'JPUTJ 

C  This  suer ou tine  o:t«':nE  total  f  !8t  tlixss  from  incut  fir* 

C  Jif. maters  ’Ey  main;  c  t  f  Tire  diameter  vs  heat  flux  t  b  1  e  (  u,  h  e  n 

u  crcvidec).  fllsc/  t  m  l  s  subroutine  is  flexible  erou;t'  to  incut 

u  tne  inn  flux  directly. 

vCM'-'CN  /IMTCH/  T  If-.'/ T  l  ,  C  ~  ,  T  1  NF 
CnfRPCTiR  *1  Pi 
ji-vf  x  -i 
PRINT  * 

PRINT  »/  "IN»L'T  -jr;  :i;.ii;tsr  (.Titers)." 

1 :  n  C  *  /  0  If  •■*  = 

PRINT  • 

irCliMF  .Li.  I  .  J )  GO  TC  IS 

L  Tne  heft  flux  CP  is  ocu  ine(  rere  trom  t  n  a  firs  d  i  ?.  "-ter 

u  vs  heat  flux  tic);-  (iihar.  crcviaed). 

PRINT  '  ( "  :ir»  Cic- ester  (meters)  =  "/E13.S)*/  SIAM15 

PRINT  *("  I, t,o  oiod  test  Hu>  (jitt5/»eiir**i!)  *  "/FI  7.5  )'/  CF 

P  k  i  N  T  * 

'unlliUS/  *)  "  C  I  A  V-  r  CF" 

v.XlTi(c5/  '(idi/ili.;))')  C- 1 A  M.  F  /  LF 
■.fiiTsCtS/  ♦) 

G C  TO  cC 
1 3  CCNTINLe 
Sri  -  *  1  ,  u 

PRINT  */  "INPUT  Imposed  meet  Flux  (u,  a  1 1 S  /  me  t  er  *  *  2 )  .  ” 

R  E  P  C  *  /  C F 1 

I F ( C F 1  , GT .  O.G)  CF  =  LF1 
FRINT  » 

PRINT  '("  Imaossd  rest  Flux  (uiJtts/ir.«t»r**i)  =  "/ST3.5)'/  C: 
PRINT  » 

V.RlTc(0ir  *)  "  ClflMF  CF" 

WRlTcUS,  '("  3X/E13.5)')  CF 

WR1T6U5/  *) 
cC  CONTlNLc 

PRINT  */  "PRESS  PtTLF.N  to  Start  Prc;ram." 

R  t  A  0  *  (  A  1  )  *  /  if 

RETURN 

eNC 

SJi  ROUTINE  1MTV4L 
PARAMETER  (  I0-=1CC  ) 

COMMON  /GECMCM/  rlL/XLL/AL/AW/Ae,V/VT 

common  /inxtck/  t:^ZxT1/.=/t:nf 

COMMON  /  T  I  y  u  M  /  TIMPXr  :Tl;<L/CT  IPG/ ST  IMPS,  TIPLrTIMGrTIMPfi 
C  CMvON  / NC v /  N 

COMMON  /  VLCM/  TWi/TLS  (ICM)  /  TCCCM),TPC(1  C.--.+  1  )  rTINFerTINFS 
CuPMuN  /WLC.M/  if 3/LKRC6# jCLS/OFS/OSaOSC OM)xSCLS(ICM) 

COMMON  / V G C M /  f  WG/T  C-/ XMV/ XMA/TlNrT 

COMMON  /WGCM/  'If  T,CRRCrSCC-rFf  xPVxF#CRPI/ .RLxXL6FF,CXM,CCL 
COMMON  /v  »CV/  :(■: :  ;;m.)  y  PST  (C  :  IGM)  ,l(G:  ;  CM)  x  ZF  (C  :  ICM  +  1  > 

CCMMON  /TO^Cy/  FSTTCP/vi-JOTCPrTCT/TI/Tc 

timl  *  tim: 

TwS  -  TZ 
00  1u  x=1/N 
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TkLCN)  =  T l 
7  C  t  n. )  i  T  I 
1l  caiil.L: 
i  >;■  =  km; 

1  'At  s  T  4 

1  u  -  TZ 
TCI  *  K 
C'XN  =  C.O 
'.ri  =  t  r 
.FT  =  wr 
Crs  =  CF 
T  IN  F  $  =  KNF 
TINFT  i  TINF 
T1N:S  *  TINF 
CALL  IMM 

i  ( iv )  =  rl  *  h  L  /  N 
2C  CCMIM.Z 
K  cTUSN 


ilrkul^:  :  ‘i :  T  V 

UAmc:*  /?m75Cm/  ?/C / -'-Tm 
C  •;  X  ’••  C  N  /  3  =  C  y  C  M  /  ML/XLL/AL/AV.  /  ii/V/VT 

V  C".,-CN  /FAT  r  •-  F  ■:  iv  /  XF  ;Lu/X''CLV  /  RFC/AriC  If  t  ;  Ti,  XLV,  XL/T1  /  PI 

C  *0N  /  v  GC  f. /  T  w  3 /  7  G  /  X V  /  <  V 4  /  T I S p  T 

C  3 1' /  V *  /WvCMX  /uC  3^Pi/pV  /  P/  ’.'kH  ,;SL  /  X^ZFF/'XW/QCL 

f  v  i  F1*:a?(  (XN0LV*XLV/.O«<1  .  C  /  7  1  -  1  .  C  /  T  G  )  ) 

xt’V  i  Cx»CIV»\//k)»FV/7(. 
a  /-l  fi  -  {(MiLi»V/x)*(F«11  -  V  V  )  /  T  Cj 
k  c  7  U  k  \ 

;  NL 

r  k  C  j  k  A  X  T  A  N  x 
Pkf.AMtTsk  C  I  j  X  =  1  C  C  ) 

CCKKON  /1MTC*/  7  JP4  ,  TZ/CF  ,7  l  Nr 

CCXX-CN  /7 IXCM/  TIi“XX,:TlML,CTIXG/G7IMPR/TlK.L/TIX3/T!MDR 
CC'XXCN  /NCX/  N 

CCXMCN  /VLC V/  Tv:/U;  CIX.KTC  <:CM),T?C(IDM+1  >/T:NF?pTINFS 
CCi'XON  /  W  L  C  x  /  ;Pi/i.SK'i/:;CLE'/CFS/«RROS(lJM)/:CLS<lCV) 
CCXMCN  /  V  G  C  FI  /  ThGrTG/XMV/XXA/TINrT 

CCXXCN  / wGCX/  .FT/wRRC/iCG/PA/PV/P/CRftlpC' RLrXL5PF/CXM/CCL 
CCXX.ON  /TOFCM/  rSTTC?/Cc!;OTlP,TC1/TI,7e 
CCXMCN  /  C  V  C-  C  M  /  CTxC/iTG/GXPVK'XPA/DXMO/CXMl 
PF'SI  =  1.45C45--* 

CALL  INFUTT 

call  input; 

CALL  INFUTI  ■ 

CALL  INITVAL 

C  r  cN  C  65x  HLs»*20CCC“.;<«0,JTL' ) 
oPZNCctp  PlLea'iOCCCnAB.JwTG') 

C  FiN ( Of /  FlL?»'2CCCCnAR*0UTl') 


aRITfc(c5p*)" 

<.  T  C  ( 1  ) 

*  FIT  :  (  6  i  /  *  ) 

t  :  L 

T  C  ( N  ) 11 

T  «  S  ( 1  ) 

TkSCN) 

wRITc  <cZ  , * <" 

t  :  j 

Tv»  G 

TG 

C  XXV 

6)  *) 

uRITctci/*) 

X  X  A 

P  (F  SI) 

C  X  F-.  I 

CXViO" 

kRITc (e7/  '  C” 

PST70R 

CrNuTC  P 

TCT 

6  Ts 

LRlTi(67,*) 
TiXPk  -  -C.CC1 
1C  CCNTINL3 

PRINT  */"  a 
CALL  WL 

TWc 

*  T  1 V  Z 

T  i  ML 

TIMG") * ) 

90 


C  tkiM  <i"  5  " 

call  Ovc 

L  PAINT  «/"  t  " 

.;'Ucs,'(5CX,:1;.i))')  T  C  P‘  L  /  TMS(1>,  TUSCN),  T  C  C 1  )  /  TCCN) 
«Kn:Ctt/*(;(2x#"1:.s))*)  Urj/  T  T  j  ,  XP,V,XML/:>*°PS:,0XPI/3XP'0 
*k:Ti(e7/'<7(0X/i1;.n)')  fSTTC.Or0iNCTCP/TCTyT9i-T«r^TT''l^TIM5 
T  i P  P  «  =  !;.V'  *  3  T  i  v *  L 

iPltlNFt  7  IKM))  STC£ 

c\.  CC-MiMi 
0  ^ifll  •/"  ?  " 

CilL  GM.T71 

tips  =  t  i « +  ct:  mc 

ifcriv'/  .it.  (f.vLH.C'.wrnHn  oc  to 

: ;  cct. 

C  P^INT  :  *• 

Call  ii.cr 

L  PRINT  -/••  r  " 

C  4  L  '■  C  r 

C  PP1KT  !  :  •' 

C-LL  lnLTT: 

c  p ■; i \i  * f "  ii  *' 

Call  RJcCC'i 
u  - 1 :  \  t  * , "  1 2  " 

KH  =  T  1  L  ♦  :  T  :  M L 

I-(TiML  .LT.  nini-C.k'Jl^Tll'-i))  C-C  To  I-., 

i.  “  (  T  I  P  L  .LT.  T  1  P  P  F, )  ,C  TO  c  0 

CO  TO  1C 
„\C 

SUcPjUTINt  sUCY 
PtSix-T  t«  C  I  3 v-  =  1  C C  ) 

C  C  P  P  C  N  /  ?  n  Y  S  C  M  /  :/i/«  A  Y  P 

C  CPKCN  /  MkT  p:  9  CP/  XPvL“/X'JCLV,khCr5pOi,£rTL,xLV/-XL/'T1,P1 

COMMON  /‘liST'LCf'/  RCCp/cFSI/Cl/CZ/tixCS/CC-ACcV 
C CP  PON  /P.CV/  N 

CLP  NON  /  V  L  C  n  /  Tv.a/T.iClLPO/.TC  (1CM)/TPC<  10M+1  >  /  TINr2/TIN*S 
CON  non  / v  Y  C  n /  :<C:ICP)/PaTlC:I0K)pL<C:IC*O,zf<C:I0P  +  1) 

C  Op  NON  /  T  0  P  C  N  /  P  S  T  T  C  P  /  If.  NQTOP/TCT/TIaT? 

P  0  T ( 0 )  5  L.C 
jC  K  <  -  1  /  N 

0  L 1  a  =  >0:T-'.*(T«o(<)-TC(K)) 

C  ILi  -  C:*C(N)-Z(P-1  )) 

p$Kk)  =  (  ? st c < - 1  >•»  to .v/4.o  +  c:o:*clta*. u  .c/3.0)  >•*<«. o/i.o 
1C  CCNTiNCi 

P  S  1  T  0  0  =  P  ;  T  (  N  ) 

uLTa  =  G»jf.  T4«(Tv»i(N)“TClli)> 

»-.notop'  =  C2*;lt«»fsttcp*»o. 0/4.0) 

A  i  T  0  K  N 
:  NO 

SJtRCOTlNc  CCNc 
PAAAMeTSA  (  IOP  =  1u  ) 

CCPNCN  /CfCPCM/  HL,XLL/tL/iV./A&/V,\<T 

C  CP  NGN  /.vlTPQPCM/  XP  vLi/  X  MC  IV/khO  /  FHOZ/  :STt/XLV,XL/T1,P1 
CCPPCN  /riATCPCP/  CP/CPLaCVP/CVV,C4V,CPC- 
c cp non  / t i y c v /  t irp x,ct :r'i, :t:pC/ct in?k» t:plaTiv",/Tikpr 
C  C  M  v  0  N  /'.CP/  !, 

C  ON  NON  /  V  L  CM  /  T  U  'i ,  T  C  >  (  I  0  M  )  ,  T  C  C  I C  P. )  /  T  PC  (  I  C  p  ♦  1  )  /  T  I  N  1  ?  /•  T  I N  e  S 
COMMON  /  C  V  L  C:  /  CT.,:  ,0TL«  (::«)  ,CTf  C 

C  CP  MOM  /„lCP./  C  r  -,C  ;  .C-  /  -.CL1  /  IrS,  C*40S<  ICO  /  OCL  $  t  I  CM  ) 

COP  MON  /  v  Y ■: m /  C  C o  :  i  C.N )  /  p  ST  c c :  i »« ) , c  (  c :  i •:  M > ,  z  p  cc : ;  jp ->i ) 

CCPMCM  /TCP  C-'/  PSTTCP/ -rNOTCP/TCT/1  i/To 
CTPC  =  (4i/(PrtC*C-»VT))-;CL? 

TCT  p  UCN)  *  "  T  T  >!  l «  C  T  P  C 

I l  P  S  T  U  P  .OS.  1  .  C  L  -  :  L  )  TCT  =  TCT+(wPNC'TOF/PSTT0F)/CP 

uu  1  u  K  -  1  t  N 
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t?cu)  *  tc<<)  ♦  ;i r?«; 

1c  CcMlNLi 

T  Pc  <N  +  1  )  *  TCT 
o  <  C  )  *  C  .  0 

.?(•:>  =  :<c> 

v  0  t  -  K  8  1  /  N 

OU)  *  -F  ST  CO«XLL/  15mC:*JL> 

:f<iO  s  :  <  * )  ♦  ct;m«u<io 
cc  CCNTINLc 

ipcnm)  -  z(n) 

j  \j  .  G  K  *  1  /  N  *  1 

lr(;?(K)  .LT.  'C(H-D)  GO  TC  iC 
i  c  C  0  i\  T  1  it  C  t 

<t  ;  1  t,  a  ;-i 

ti  C  wi.  T  I  itL  t 

*^:T:(t5#«>  "$*CF  IK  5  U  -  -  '.  L  T  I  M  *  c.-il.  *»ClO  NCI  IN  CP  e  4  S I NG  .  " 
il-.s 


i  C*  -  C  C I  I  i  ^  “  t  C  C  * 

P-f.-y.iTti  (  :ov=1CC  ) 

/ > . ; n / 

c me-/  Tw-/T»ic;i-)/T(  cil v),tpccom+i >,t:n??/Tin'' 

l  - ■'*  C K  UVZt'f  I  (1!:  :ZM),KiT  (t  :  I-?')  /  L  (0:  ICH)/ ZP(C  :ICV  +  1 ) 
v  ~  w 
I  •’  8  C 
K  cc'NTlhli 
\  -  K»1 
I  G  -  I  f 
iC  COMING: 

I  f  (In  •  G  c  .  ft )  j  C  T  C  : 

Ir(iftlfM)  .GT.  Uf.))  CC  TC  ;C 
I-  =  I c  + 1 
c  0  TO  i C 
00  CCNTlNLr 

If(In  .Ne.  IG)  GO  TC  41. 

TC<fQ  -  TPC(Iri  +  1) 
o  C  TC  >  C 
4  c  CONTI  It  Cc. 
a  =  0  .  C 

J  ~  *  u 1  c 

5L  COMING* 

Ir(J  .GT.  Ib)  GC  TC  sC 
a  «  a  ♦  (I"  (J)*I°C-  -1 )) -T?C  (J) 

J  -  J*1 
CO  TO  SC 
cC  CONTINUE 

a  8  a  ♦  (Z  =  CG+1 )-:  (K-1  >)*TPCCIS-M>  ♦  ( l  U)  -  IP  (IH)  )  *T  PC  U«  +  1  ) 

i  c  c  it  j  -  j/ ceu>-:<k-i ) ) 

r«  C  ON  T I N  L  E 

IrU  .LT.  N )  GC  TO  1C 
4  t  T  c  k  N 
c  nC 

atsfouTlN.  hutta 
PfiriAMtT  £  R  (  lObsICC  ) 

Cob  MON  /TIVCM/  TI«MX/I.TI*L/CTIMC-/0TIN,.PN/7lML/'T:f'G/TlM?R 

CCNMGN  /  NO v /  N 

COMMON  /  V  L  C  v  /  VUUMI) 

COMMON  /GvLCM/  CVUCKd) 

ClNiNalCN  V1U?i«  +  1),  SCI3V.+  1) 

M  =  NM 
b  *  CTIMl 
Cl.  10  N  =  1  - v 
vl  U)  -  vco 
1c  CCMIUCr 
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(.ALL  OVL 

00  20  n  *  1  ,  v 
$oo  -  :v(K) 

V (  N)  *  VI <K)  +  CV(K)*n/2.0 
2C  CCNTINLc 
Call  Ovl 
0 C  30 

iu)  «  sen  ♦  2.0*cvC(O 

vco  =  vi  to  cvu>*h/;.g 

30  CGNTlNLi 
CALL  OVl 

cc  ;o  k*i ,y 

SU)  =  SU)  ♦  i  ,C»lv(iO 
VIA)  *  VICK)  ♦  C  V  (  K  )  «  H 
4C  CONTlNLt 
CALL  OVL 
CO  3C  K*1,M 

VU)  S  V  1  (  K  )  ♦  (S(K)  +  ;V(lO)*P/c  .c 

H  CCNTlNLc 
A  *  T  0  k  N 
e  NC 

JlicAOuT  i  N  E  OVL 
p  a  k  a  m  S  T  tk  (  iom=1CC  ) 

COMMON  /hE«T?LCM/  kCDn/:P$I/C2/C2rCB/CS»CG#C£V 
COMMON  /NIC*/  N 

COMMON  /CVLCM/  CTwt/OTCS (ICM),C1PC 

COM.MCM  /  W  L  C  v  /  IPS,  tRR0s/iCL6/Lf  S /ORfiOS  ( ICK)  ,CCL  S  (IC« ) 

Call  WL 

OTh;  =  („r?  -  'RROE  -  0  C  L  i  )  /  A  C  u» 
u  0  10  K  =  1/N 

0  T ».  S  C  K  )  A  (  C  F  5  -  CKKOSU)  -  ;CL$<K>  )/feCDL 
1C  CONTINLc 
Kcl L«N 
-NC 

MjfckCLTIN:  *  L 
PARAMOUR  (  I  D  M  *  1 C  C  ) 

COMMON  /h£«tflCm/  «CCL,E?$  I>CI/C2,C5/CS/C0/CfcV 
COMMON  /NC*/  N 

COMMON  /  V  LC  M  /  TW9#1US1ICM)/TC<10M)rTPCCI0M  +  1)»TINF8/TINFS 
COMMON  CF5/CRkC3/CCLi/CFS/vRS0S(IDM)/CCLSaC«> 

.kkOs  *  £PSI*(Tv,3**4  -  TINFE**4) 

TOlFr  *  MAXCO.C/  Twe-TC(D) 
wC  ls  =  C5*TC!PF**(4.0/3.G) 

OC  10  K*1,N 

CRKC3  (K  )  -  =PSI*  (TVSClO **4  -  TInF£**4) 

IClrr  s  MAXCO.C*  1 k S < k > -1 C <k > > 

LCCSOO  =  C$*TOIFF**U.C/i.C> 

1„  CONTINUE 
k  c  T  0  R  N 
ENC 


iUikOUlIN;  3XCITT4 
COMMON  /V3CM/  V (4 ) 

COMMON  /CVC-CM/  CVU) 

COMMON  / T I M  C  M /  TlMMX/OTlMLrCTIMG^OllMPfc/TIML^TIMG/IIMOR 
OlMiNSICN  VI (4)/  3(4) 

h  s  OTIMO 
OC  10  K  =  1  / 4 
VICO  *  VOO 
1C  CCNTlNCt 
CALL  OVC 
00  2U  A*1,4 
SU>  *  C  V  C  < ) 

V  (X)  =  V  1  (  M  )  +  CVOO«r*/< 
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t C  CONT IPi : 

Call  c vc- 

uC  Jb  X  3  1  /  4 

i  <x )  *  i  ( s. )  *  i.C*[V(M 

v<n>  •  v i ( k )  *  :vu) »)-/;. j 
:c  coMifai 
CALL  ovc 
C  0  4  U  X  3  1  /  4 

su)  *  $<o  *  :.o*:vck) 

V  (O  3  vl  (O  *  CV(K)*f* 

* C  CCxT.I.  Lc 
CALL  OVO 
L  u  5  U  K  *  1  »  4 

v c k )  *  vico  ♦  (s<K)  +  :vc'))*H/t.c 
>C  C  Of.  T  i  N  L  C 
k  c  1  U  k  N 
;NC 

SUtSCUTI.\c  C  VC 

LCPFCN  /GLCvCM/  *l/XLL* A  L  /  A  k  /  A  S  ,  V  ,  V  T 
C'.PFCN  /r-.iiTCPCP/  L>-^CkL/CVP,CVV#CA  V/CFG 
COPKCN  /PeA7slCM  kC**/iPiIyC23C2>Cv#CS/CG,Cev 
cckkoi*  /vcc*/  ThC/ic* xmv^xpa^t inpt 

CCPMCN  /CVC-CV/  CT„G*‘jTi,GXJ'V»OXMi/CXMC/CX«: 

LONDON  / wGCH/  CAT/Ck^C^CCG/PA/PV^P/CkSl/CkL/XLerF/CXM/CCL 

CCM-O.N  /TO*C*/  P$nCP/v;  NOT C ? , TC T / T  I , T r 

CCFMGN  /  P  H  Y  S  C  M  /  C/G/PATF 

C CPMG!\  /ViMCM/  AVM^CCV'iT/AVN'T/GPVNT 

CALL  W  u 

GTViG  =  (JFT  -  CPRC  -  .PkX  -  CCCO/SCCW 

jxmsn  -  al«:xk«(cvp»t:  -  cvv*tg) 

GTG  3  (CCi«AW  -  ;CL»AL  *  LXPTFV)/<XPA*CAV  ♦  XMV*CVV) 

GXP'  3  AL-CXM 
C  xp  C  3  i.  •  0 
c'xpv  s  cx«: 
l  x  p  a  ~  C«w 

Ir(?  ,Lc.  FV-.T)  «  1 1  L  A  N 

uXPC  3  A  VM  *CC'VkT«Si;k  T(  i  .  > C P-PATP) * ( XMV*XPA ) /V) 

I  r  (  P  .LT.  <?VNT  +  0  =  VXT  ) )  CXPC  *  C  X  MC*  <  P-r  VNT)  /  OP  VN  T 
CXPV  3  CXMV  -  PXffOXPV/ (XMV+XMA) 

*;  XP  A  3  -CXNC‘XMA/ 

lTo  =  C  T  G  -  ;xMC*P«V/(  (  <VV  +  XM4)*<XM»*CAV*XMV*CVV>  ) 

*  C  T  U  kN 
cNC 

susftOLT  CNi  vi Gr 

CC.V«GN  / ?iY  5  CP/  A,o,AAr.‘- 

CCP-PCN  /GcCHC*/  -L/XLL/Al/AX/A'/V^VT 

COP. VON  /-’AT ?  F  P  C  P  /  >PCLA,XP.CLV/AFO/k:-OZ/o:TS,XLV/XL/T1,P1 
CCPFON  /ntATCPCM/  CP/CPL#CVP,CW,CtV,CPC- 
C  CP  PON  /neiTPLCP/  ACLx/£PGi/C2sC<sC$*CS*CGrC£V 
CCPPON  /VGCM/  Twi/ltt/ XPV/Xf A^TINPT 

CCPMON  /  W  5  C  P.  /  CPT/LSkO/'.Ci^PA/rV/P/CftSI/CRL/XLiFS^CXK/CCL 
CCPP-ON  /7CPCM/  ‘>$TTGr',C;:NOTCP,TCT#1  !,Tr 
«|khC  3  c?SI*(TWG*«*i  *  TiNrT**4) 

TOIFP  3  P  A  X  <  0  .  C  *  T»C— To) 
wLG  3  CG-TCIFF**Ci.C/ i.C) 

PA  =  (XPA/XPOLA) »(i/V)*Tt 
rV  =  (XPv/Xy.CLV)*<R/V)-Ti 
r  3  PA  ♦  p  V 

To  3  1.C/C  ( 1  .C/T1 )  -  (r/ (XPGLV*XLV))»lCv,C?/P1)  ) 

call  GLSPT'.yP 

wRRI  3  c?5I*(TwC**a  -  T>*4) 

CkL  3  CkRI*AW/il 

XLcrF  3  XL  ♦  CPLM7I-KT) 

CALL  :VAP 
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.CL  *  C  *  .X  *  X  :  F  c  -  (  .  •,  L 

*  c  T  L  K  N 

cnC 

iUC*OL)TlN;  iVJP 

CC KMCN  /HtaTCPC"/  Ce/CPL»CV»,Cvv/CiV/CPu 
LCNMON  /'■£ATsLCf’/  FCCh/rBSI/CI/C</C?/Ci»CG/CEV 
COMMON  / TCL  C'V  T0F1CL/XTCL/FTCL/NTCL 
CCi'VCN  /VoC"/  T»G/1C,XMV/Xf'4*TlKF  T 

CC?  f*CN  /WCC'"/  CFT,LFSCxCC5#Pfl/PV/P/CRRI/C|KL/XLtFF,CXN/CCL 
CGFi-GN  /TO'Cv/  PSnCP/CFNOTCS/TCT/TI/Te 
T  j  I  F  F  =  K  i  X  ( 0 . 0  ,  T  C  -  T  I ) 

a  *  cev»Tc:FF».n.c/i.o> 

®  *  COC-*TOICF/XL=FF 
C  =  .  R  L  /  X  U  r  F 

:fct;:ff  ,c-t.  tcftclj  gc  tc  io 
i.ip  =  c 

A  C  1  0  n  N 

1C  v-CNTI/ViCc 

Call  £  V  I.  E  w  T  <  J  ,  E  ,  O  C  x  r ,  1 1  F  ■_  »  I  ?  6  L ) 

:  F  C  <  i  1  F  L  ♦  ;  r  u  .LT.  NTCL)  RE TJRN 

h*lT;(tJ/  *)  "il'P  IN  it.ACCTINc  EVAP.  NT C L  EXCEEDED." 

i TC  F 

e.NC 

ilie'CjT  IN:  “  V  Nix  T  (4/  F  /  C  t  X/  1 1  F  L  /  I  £  F  L  ) 

CCr'WuN  /TCLCM/  TOFTCL/XTCL/FTOL/NTCL 
IF  (A  .  GT.  xTOL)  T  h  fc  N 
AVAR  =  A*SCRT(5*C/ (A»A*c  *i*C)  ) 
x  s  c  +  max(avar/C.fv*ft;l) 

Cl.'C 

X  *  C  ♦  o.Si*?TCL 
i >; C  IF 
I1H  =  C 
lirL  =  C 
;c  £0/  n*1/ntcl 

r  *  X  -  A  *L  C  j ( 1  .  0*  t  *  X  /  ( X-C )  ) 

of  =  i.c  ♦  a*(  <  s*c/((i.c+s)*cx-c)+e*c)  )/(x-c>  ) 

CX  *  r/CF 

IF  (F  .GT.  C.O)  TrSN 
1 1  F  L  =  1 1  F  L  ♦  1 

IF((x-C)  .Cl.  FTOU  x:TLKN 
CX  =  M  1  N  (  OX  /  0.  5*0-0) 
cLSc 

I2FL  -  12  cl  + 1 
ENC  IF 

IF  (  A  3  j  (  r  )  .IT.  -TCI)  4-TUH. 

X  =  X-C  X 

Co  ccntinle 
k  E  1  CRN 
E,\lC 

iciSCUT  IN-.  <CR?T:fF 

v.CF'ICN  /VGC  lM  T  W  .4  /  1  C  <  X  '<  V  /  X  F  A  ,  T  I  N  F  T 

CCN>'On  /XGCW  CFT,CACC,CCO/PA,rV/P/CRLOCRL/.XLEFF,CX«,CCL 
COMCN  /  T  0  c  C  v  /  FSnCF/.£.'.CTCF,TCT^TI,T-: 

t;  -  n  n  ( o .  5  ♦  (  t  c-  -  t  c  t  )  /  t  •: ) 

A  t  T  Cf  k  N 
c  i'.L 
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APPENDIX  C 
EXPERIMENTAL  EFFORT 


Specially  designed  fire/tank  tests  were  conducted  to  measure  the 
response  of  a  tank  containing  liquid  to  fire  impingement,  ihis  experimental 
effort  consisted  of  measuring  effects  of  the  heat  input  from  a  large-scale 
turbulent  sooty  fire  on  a  partially  filled  liquid  solvent  container.  The 
goals  of  this  task  were:  (1)  to  evaluate  the  fire/tank  interaction 
environment,  and  (2)  to  obtain  the  data  necessary  for  preliminary  validation 
of  the  numerical  computer  prediction  model. 

In  these  tests,  commercially  available  55-gallon,  18-gage  steel  drums 
were  used.  The  drums  were  instrumented  to  monitor  the  internal  pressure  and 
temperature  rise  of  a  fire-engulfed  tank.  The  drums  were  equipped  with  a 
network  of  chromel  alumel  thermocouples  for  measuring  the  variation  of  the 
liquid  temperature  (Figure  C-l).  A  pressure  gage  (0-30  lb/in2)  was  used  to 
monitor  the  pressure  rise  resulting  from  expansion  of  the  vapor  (Figure  C- 
2)  and  manually  vent  the  tanks  before  rupture.  The  tanks  were  placed  about 
0.15  m  above  the  surface  of  a  4.3-m-diameter  JP-4  fuel  pool  fire  in  both 
vertical  and  horizontal  orientations.  Prior  to  the  start  of  the  fire,  the 
tanks  were  filled  with  water  up  to  12  centimeters  (5  inches)  from  the  top  of 
the  tank.  Hie  temperatures  were  read  by  a  data.ogger  at  5-second  intervals, 
and  each  event  was  recorded  conti nously,  using  video  and  still  photographs. 
Figures  C-3  and  C-4  show  the  fire  pit  and  the  orientations  of  the  tank  witn 
respect  to  the  pool  fire.  Figures  C-5  and  C-6  show  the  temperature  profiles 
for  two  tests  with  the  tanks  oriented  vertically. 

General  analysis  of  the  temperature  data  and  postfire  inspection  of  the 
tank  shell  revealed  the  following: 

1.  The  tank  initial  pressure  rise  Is  predominantly  dua  to  the  expan¬ 
sion  of  vapor  in  the  ullage  volume, 

2.  The  radiation  heat  emitting  from  the  dry  wall  in  the  ullage  vo; .me 
is  the  main  heat  flux  mechanism  in  the  gaseous  phase. 
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from  drum 


Figure  C-2.  Pressure  Monitoring  Assembly 
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(a)  Vertical  orientation  tank  fire 


(b)  Horizontal  orientation  tank  fire 


Figure  C-3.  Tank  Fires  with  Different  Orientations 
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(b)  Horizontal  orientation  test 


Figure  C-4.  Vertical  and  Horizontal  Orientation  Tests 
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Tank  Temperatures  (°F) 


Figure  C-5.  Fire/Tank  Test  No.  1  Temperature  Profiles 


Tank  Temperatures  (°F) 


Thermocouple  tree 


Figure  C-6.  Fire/Tank  Test  No.  2  Temperature  Profiles 


Fire  Temperature  (°F) 


3.  Discoloration  of  the  wall  next  to  the  vapor  compared  to  the  wall 
in  contact  with  the  liquid  is  indicative  of  a  large  temperature 
difference  across  the  vapor/liquid  interface,  and 

4.  The  thermocouples  placed  at  various  elevations  along  the  height  of 
the  tank  showed  a  decreasing  temperature  gradient  downward  from 
the  top  of  the  tank. 

The  discoloration  of  the  tank  shell  suggests  a  strong  conduction  heat 
transfer  in  the  tank  wall  at  the  liquid  meniscus  which  contributes  to  the 
overall  thermal  activity  of  the  interface  and  should  be  considered  in  the 
model . 

The  predictions  of  the  model  proved  to  be  consistent  with  physical 
intuition  and  agreed  well  with  the  overall  trend  of  the  experimental  data. 
However,  the  predictions  from  the  calculations  consistently  showed  higher 
values  than  the  corresponding  experimental  data.  This  variation,  sometimes 
as  high  as  200  percent,  persisted  during  all  test  events. 

The  factors  that  might  cause  or  contribute  to  the  discrepancy  between 
the  simulation  predictions  and  the  tests  were  evaluated.  The  assumption  of 
constant  heat  flux  over  the  entire  tank  in  the  simulation  may  not  have  been 
duplicated  in  the  tests  as  a  result  of  wind  effects  and  the  overall  test 
configuration.  Also,  because  of  safety  considerations  the  tank  was  vented 
manually,  which  could  have  affected  the  tank  response.  Project  resource 
limitations  did  not  permit  further  testing  and  validation  of  the  tank 
response  model . 
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